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Abstract 

A  Comprehensive  Study  of  Vortex  Breakdown  Flow  Mechanisms: 

Computational  Investigation  and  Preliminary  Control 

by  William  F.  Johnson  III 
Chairperson  of  the  Supervisory  Committee 

Professor  Mitsuru  Kurosaka 
Department  of  Aeronautics  and  Astronautics 

A  modified  vortex  filament  method  is  used  to  simulate  the  evolution  of  the 
transient  formation  of  vortex  breakdown.  The  method  supports  previous  studies, 
illustrating  that  vortex  breakdown  is  initiated  by  a  negative  vorticity  gradient  which 
triggers  an  inviscid  self-induction  feedback  mechanism  and  when  subsequently  subjected 
to  viscous  effects,  results  in  steady  state  vortex  breakdown.  The  results  of  the  method  are 
first  validated  experimentally  with  numerous  past  dye  flow  visualization  and  particle 
image  velocimetry  investigations,  and  then  used  to  qualitatively  investigate  the  self- 
induction  flow  mechanisms  during  the  formative  stages  of  transient  breakdown.  As  a 
complement  to  the  qualitative  investigation,  a  quantitative  analysis  is  performed,  which 
yields  a  local  and  dynamical  relationship  relating  the  azimuthal  vorticity  gradient  at  a 
particular  location  to  the  curvature  of  the  instantaneous  streamline,  projected  onto  the 
meridional  plane,  at  the  same  location.  This  relationship  further  shows  that  once  radial 
expansion  commences  in  the  region  of  negative  azimuthal  vorticity,  it  continues  to 
expand  such  that  the  meridional  streamline  becomes  more  curved  with  time,  supporting 
that  the  negative  vorticity  gradient  not  only  initiates  the  radial  expansion,  but  also,  feeds 
its  subsequent  growth.  On  the  contrary,  in  the  region  of  a  positive  gradient,  the 
streamline  continues  to  flatten  fostering  radial  contraction  of  the  vortex  tube,  which 
provides  a  closure  to  expansion.  In  attempt  to  suppress  breakdown  in  two  preliminary 
control  simulations,  this  positive  azimuthal  vorticity  gradient  is  then  introduced  to  the 


vortex  flow  just  prior  to  breakdown.  Results  from  these  control  simulations  illustrated  a 
temporal  and  spatial  delay  in  breakdown  as  well  as  exhibiting  flow  behavior  associated 
with  complete  elimination  of  breakdown. 
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1  Introduction 

1.1  Introduction  to  Vortex  Flow 

Vortex  flows  are  regions  of  swirling  fluid  which  can  occur  in  any  type  of 
medium.  Examples  of  vortex  flows  are  apparent  in  our  everyday  lives,  from  complex 
weather  phenomena,  such  as  violent  tornadoes  to  the  simple  swirling  water  pattern 
formed  when  draining  a  sink.  Whether  one  is  driving  behind  a  tractor  trailer  or  piloting  a 
high  performance  aircraft,  one  can  not  escape  the  presence  of  these  vortex  flows. 

In  the  area  of  aerodynamics,  vortex  flows  are  also  commonly  found  emanating 
from  and  around  aircraft,  specifically  aircraft  marked  by  large  wing  sweep  and  leading 
edge  extensions.  Figure  1-1  illustrates  an  F/A-18  aircraft  with  a  leading  edge  extension 
generated  vortex  while  Figure  1-2  illustrates  the  structure  of  a  similar  wing- tip  vortex. 


Figure  1-1:  Vortex  flows  over  F/A-18  {ref.  1) 

Leading  edge  vortices  such  as  the  one  illustrated  over  the  F/A-18  are  sources  of  high- 
energy,  relatively  high  vorticity  flow  which  generate  a  strong  component  of  axial  velocity 
along  the  core  of  the  vortex.  The  local  static  pressure,  which  in  general  becomes  lower 
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toward  the  vortex  core,  is  further  reduced  in  the  vicinity  of  the  strong  axial  velocity  and 
near  the  leading  edges. 


Figure  1-2:  Wingtip  generated  vortex  (ref.  2) 


This  reduction  in  pressure  creates  a  strong  suction  on  the  top  surface  near  the 
leading  edges,  enhancing  the  lift  while  energizing  the  flow  in  the  boundary  layer  to  help 
delay  separation.  For  this  reason,  the  lift  coefficient  curve  for  a  delta  wing  exhibits  an 
increase  in  lift  coefficient  for  angles  of  attack  at  which  conventional  wing  configurations 
would  be  stalled  (ref.  3).  Unfortunately,  this  increase  in  lift  performance  comes  at  the 
price  of  increased  drag  and  required  thrust. 

At  very  large  angles  of  attack,  however,  these  vortex  flows  experience  a  dramatic 
form  of  flow  disruption  referred  to  as  vortex  breakdown.  This  phenomenon  is 
characterized  by  reversal  of  axial  flow  and  the  sudden  explosion  of  the  vortex  core  (ref. 
4).  The  breakdown  poses  severe  limitations  on  the  performance  of  an  aircraft  due  to  the 
sudden  degradation  in  aerodynamic  forces  and  moments.  Moreover,  the  rapid 
fluctuations  within  the  turbulent  breakdown  region  can  cause  buffeting  and  structural 
damage  through  high-cycle  fatigue,  especially  on  the  empennage  structure  located 
immediately  downstream  of  the  breakdown  region.  This  effect  is  prevalent  in  many 
aircraft,  such  as  in  the  case  of  the  F/A-18,  seen  in  Figure  1-1. 
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Over  the  past  few  decades,  vortex  breakdown  has  been  a  subject  of  rigorous 
theoretical  and  experimental  study.  During  this  period,  a  great  deal  of  effort  has  been 
spent  in  the  hopes  of  understanding  the  flow  mechanisms  of  vortex  flow  and  predicting 
its  occurrence.  These  efforts  have  included  an  assortment  of  experimental  investigation 
methods  such  as  flow  visualization,  small-scale  model  tests  in  wind  tunnels,  as  well  as 
full-scale  aircraft  flight  tests.  In  addition  to  the  vast  experimental  work,  a  great  deal  of 
numerical  techniques  has  also  been  employed  which  include  an  array  of  grid  dependent 
Euler-base  numerical  methods  and  more  comprehensive  Navier-Stokes  code  analysis. 
Despite  some  progress  in  predicting  the  occurrence  of  vortex  breakdown,  the 
aforementioned  studies,  along  with  current  initiatives,  continue  to  fail  to  grasp  the 
underlying  flow  physics  of  the  swirling  flow  As  a  result,  a  universal  and  general 
understanding  into  the  complex  flow  mechanisms  that  affect  vortex  flow  and  the  physical 
process  governing  its  breakdown  still  has  not  been  reached  {ref.  5). 

The  purpose  of  this  research  is  to  investigate  the  flow  mechanisms  which  occur  in 
the  formative  stages  of  vortex  breakdown  with  the  goal  of  eventually  using  these 
associated  flow  behaviors  to  employ  methods  to  suppress,  delay,  or  weaken  vortex 
breakdown.  This  is  the  first  and  most  crucial  step  in  the  over-all  effort  to  eventually 
suppress  the  negative  effects  of  vortex  breakdown  while  maintaining  and  enhancing  the 
positive  attributes  associated  with  vortex  flow. 

1.2  Vortex  Methods 

An  alternative  to  the  grid-dependent  Euler  and  Navier-Stokes  numerical  methods 
is  vortex  methods  which  simulate  the  fluid  flow  under  an  assumption  of  non-linear 
dynamics  of  vorticity  {ref.  5).  In  this  method,  the  inviscid  motion  of  vorticity  is 
determined  by  the  local  fluid  velocity,  which  in  turn  is  kinematically  calculated  from  the 
prescribed  vorticity  field.  In  addition  to  the  robust  and  simple  nature  of  vortex  methods, 
another  strong  benefit  of  such  methods  is  the  ability  to  simulate  the  real  time  dependent 
build-up  of  the  vortex  flow.  Unlike  the  many  other  numerical  methods,  the  vortex 
method  is  able  to  model  the  flow  from  the  very  beginning  and  further  track  the  continual 
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development  of  the  swirling  flow  in  its  formative  stages,  without  requiring  the  initial 
introduction  of  prescribed  vorticity. 

The  most  promising  and  robust  vortex  method  is  the  vortex-filament  method. 

One  of  the  first  attempts  at  utilizing  the  vortex- filament  method  to  model  the 
phenomenon  of  vortex  breakdown  was  carried  out  by  Nakamura,  Leonard,  &  Spalart 
( ref.  6).  In  the  current  investigation  a  modified  version  of  the  vortex-filament  method 
initially  developed  by  Cain  ( ref.  7)  will  be  used  to  simulate  the  complex  vortex  flow. 

Cain’s  vortex  filament  method  allows  one  to  look  at  the  temporal  evolution  of  the 
transient  formation  of  vortex  breakdown  with  the  use  of  the  Biot-Savart  law  to  calculate 
the  motion  of  fluid.  The  Biot-Savart  law  is  based  on  the  induced  velocities  associated 
with  a  specified  vorticity  distribution.  Therefore,  if  the  law  is  modified  such  that  it 
includes  only  the  induced  velocity  profile  of  the  vortex  core  and  the  effects  of  vortex 
stretching,  we  can  limit  the  governing  mechanisms  of  the  vortex  flow  to  those  parameters 
associated  specifically  with  induced  velocities,  the  vortex  core  velocity  profile,  and 
vortex  stretching.  With  this  in  mind,  let  us  now  explain  how  the  filament  method  is 
employed  to  model  a  vortex  tube. 

1.3  Cain’s  Vortex  Filament  Method 

Although  modifications  will  be  made  to  Cain’s  original  vortex  filament  method, 
the  discretization  and  initialization  process  of  the  modeled  vortex  tube  will  remain 
unchanged.  Therefore,  in  order  to  aid  in  the  discussion  in  the  following  chapters,  a  quick 
overview  of  the  vortex  filament  method  is  necessary.  For  a  complete  and  detailed 
explanation  of  the  method,  please  refer  to  Cain  ( ref.  7) 

1.3.1  Brief  method  overview 

A  vortex  tube  is  composed  of  many  vortex  lines.  These  vortex  lines  are  always 
aligned  with  the  vorticity  vector  associated  with  the  fluid  motion,  as  seen  below  in  Figure 
1-3.  In  essence  these  vortex  lines  are  nothing  more  than  the  simple  vortex  filaments 
which  will  be  used  in  the  simulation. 
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Figure  1-3:  Vortex  tube  ( ref.  7) 

Using  the  Biot-Savart  law,  one  can  calculate  the  velocity  induced  by  the  vorticity 
of  one  vortex  line  on  another  nearby  vortex  line.  In  doing  such,  we  can  determine  the 
movement  of  the  vortex  lines  relative  to  each  other.  Helmhotz  vortex  theorem  states  that 
the  vorticity  associated  with  a  fluid  will  move  with  the  fluid  line.  Therefore,  as  induced 
velocities  cause  movement  in  the  fluid  elements,  they  likewise  cause  movement  in  the 
associated  vorticity  elements.  With  this  idea  in  mind,  we  can  calculate  the  induced 
velocities  of  the  vortex  lines  at  a  given  time,  multiply  those  velocities  by  the  time  step, 
and  ultimately  calculate  the  movement  of  the  vortex  lines  relative  to  each  other. 

Repeating  this  process  at  every  time  step  offers  us  the  ability  to  track  the  temporal 
evolution  of  the  vortex  lines.  Therefore,  by  representing  a  real  vortex  tube  with  a  group 
of  vortex  lines,  prescribed  with  a  given  vorticity  distribution,  the  vortex  filament  method 
can  model  the  motion  of  a  real  vortex. 

1.3.2  Discretizing  the  vortex  tube 

To  use  the  vortex  filament  method,  we  need  to  discretize  a  vortex  tube  into  an 
appropriate  number  of  vortex  filaments.  In  order  to  appropriately  discretize  the  simulated 
vortex  tube,  the  vortex  lines  are  prescribed  such  that  their  sum  occupy  the  same  cross 
sectional  area  as  the  vortex  tube.  Each  vortex  line,  at  a  given  vortex  line  radius,  is 
assigned  the  same  Gaussian  distribution  of  vorticity,  similar  to  that  of  the  overall  vorticity 
value  of  the  tube.  These  lines  are  then  packed  densely  enough  to  assure  that  there  is  no 
empty  space  inside  the  discretized  tube.  Additionally,  in  order  to  maintain  a  smooth 


6 

distribution  of  vorticity  within  the  tube,  the  cross-section  of  vortex  lines  overlaps.  A 
sketch  of  the  discretized  vortex  tube  is  illustrated  in  Figure  1-4. 


Figure  1-4:  Vortex  tube  discretization  (ref.  7) 


Now  that  the  vortex  tube  is  discretized  into  many  overlapping  vortex  lines,  it  is 
further  discretized  individually  along  those  lines.  To  do  this,  the  vortex  lines  are  broken 
up  into  smaller  vortex  filament  segments,  illustrated  below  in  Figure  1-5. 


Figure  1-5:  Vortex  filament  segments  (ref.  7) 


The  endpoints  of  each  vortex  filament  segment  are  referred  to  as  nodes  and  mark  the 
location  at  which  the  induced  velocity  associated  with  each  respective  filament  is 
calculated.  In  short,  the  movement  of  the  filaments  is  modeled  by  tracking  the  temporal 
evolution  of  nodal  position. 

The  simulated  tube  is  discretized  into  nine  concentric  cylinders.  These  cylinders 
have  radii  of  0.42,  0.36,  0.30,  0.24,  0.18,  0.12,  0.08,  0.06  and  0.03.  Each  cylinder  is  then 
further  discretized  into  sixteen  vortex  filaments,  which  are  even  spaced  along  the 
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circumference  of  the  cylinder.  Additionally,  a  tracking  fdament  is  also  specified  along 
the  centerline  of  the  vortex  tube,  r=0,  which  contains  no  vorticity. 


1.3.3  Initializing  the  vortex  tube 

Once  the  vortex  tubes  are  specified  in  the  radial  direction,  the  next  step  is  to 
initialize  them  in  the  axial  direction.  In  this  step,  the  starting  axial  location  of  all  vortex 
filament  segments  is  defined  by  specifying  the  endpoint  node  locations  of  each  filament 
segment.  For  this  investigation,  the  initial  axial  filament  length  of  0.1  was  chosen. 
Therefore,  at  a  given  radius,  all  the  nodes  are  equally  spaced  axially  along  the  vortex 
centerline  at  a  value  of  dz=0.1. 

After  the  initial  location  of  the  vortex  filaments  has  been  specified  and  given  a 
core  radius,  vorticity  is  then  assigned  to  each  vortex  filament  segment.  Each  individual 
vortex  filament  segment  is  assigned  a  vorticity  value  by  prescribing  to  each  segment  a 
value  of  circulation,  TUj ,  such  that  the  total  circulation  of  all  of  the  concentric  cylinders 

at  an  axial  location  is  equal  to  the  total  circulation  value  of  the  vortex  tube,  at  that  same 
axial  location. 

Since  vorticity  remains  constant  along  a  vortex  line,  each  vortex  filament  along  a 
specified  radius  will  have  the  same  value  of  vorticity.  On  the  contrary,  the  vorticity 
distribution  among  differing  radii  will  be  assigned  such  to  match  that  of  an  ideal  vortex 
created  by  unbounded  swirling  flow.  The  circumferential  velocity  and  the  axial  vorticity 
of  this  ideal  flow  are  defined  by  Equations  (1-1)  and  (1-2)  respectively. 
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To  match  experimental  data,  Cain  choose  a  value  of  b  =  10.9.  In  the  above  equations  C  is 
defined  as  the  total  circulation  value  at  a  given  axial  location  and  is  the  sum  of  the 
circulation,  T;  . .  Figure  1-6  shows  the  circulation  in  each  concentric  ring  of  filaments  as 

well  as  the  total  circulation  of  the  discretized  vortex  tube. 
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Total  Vortex  Circulation:  Ideal  Verses  Discretized  Case 


Figure  1-6:  Ideal  vs.  discretized  circulation 

The  series  of  small  amplitude  colored  lines  illustrates  the  sum  of  the  circulation  from  the 
sixteen  vortex  filaments  in  each  of  the  nine  concentric  rings.  The  sum  of  circulation  of 
all  vortex  fdaments  (all  nine  rings)  within  the  vortex  tube  is  represented  by  the  dashed 
black  line,  and  is  further  compared  to  the  total  circulation  distribution  of  the  ideal  case, 
shown  by  the  dashed  red  line.  From  the  data  offered  in  Figure  1-6,  a  Gaussian 
distribution  can  be  obtained  to  approximate  the  distribution  of  vorticity  in  each  vortex 
fdament  core.  This  approximation  is  denoted  by  Equation  (1-3). 

2ao/r/  (rj  +  r2  iacjj) 


c-co7  =cr\ 


i(r;  +  r)2  +  aa2j  J(r;.  -  rf 


+  ao 


,  for  j  =  1  to  9 


(1-3) 


In  Equation  (1-3)  the  j  subscripted  values  refer  to  the  value  of  the  quantities  for  an  entire 
concentric  ring.  To  best  match  a  Gaussian  distribution,  we  choose  a  =  0.413.  More 
specifically,  C  •  T  is  the  circulation  in  the  entire  jth  ring  of  vortex  filaments,  and  since 
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each  concentric  ring  is  made  up  of  sixteen  filaments,  the  circulation  of  an  individual 

271-C-F. 

filament  in  that  ring  is - . 

16 

1.3.4  Prescribed  Vorticity  Gradient 

In  his  investigation  of  transient  vortex  breakdown,  Cain  concluded  that 
vortex  breakdown  was  brought  about  through  a  self-induction  caused  by  a  negative 
vorticity  gradient  in  the  vortex  core.  With  the  use  of  his  vortex  filament  method,  he 
showed  that  this  negative  vorticity  gradient,  with  azimuthal  vorticity  in  the  vortex  tube 
increasing  in  the  upstream  direction,  was  an  important  initiator  of  the  vortex  breakdown 
mechanism.  For  this  reason,  we  also  want  to  highlight  this  negative  vorticity  gradient  in 
our  model  and  will  do  so  by  causing  a  change  in  circulation  in  the  vortex  tube.  In  order 
to  create  this  negative  vorticity  gradient  case  (vorticity  in  the  vortex  tube  increasing 
upstream),  the  circulation  profile  in  the  vortex  tube  coming  out  of  the  vorticity  source  as 
a  function  of  time  will  be  linearly  increased  from  an  initial  and  minimum  value  to  a  final 
and  maximum  value. 

1.3.5  General  Algorithm 

Figure  1-7  shows  a  brief  flow  chart  outlining  the  general  algorithm  of  the  vortex 
filament  model.  The  complete  and  modified  simulation  code,  written  in  C++  is  located 
in  Appendix  D. 


Figure  1-7:  General  algorithm  (ref.  7) 
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2  Preliminary  Modifications 

The  self-induction  theory  of  vortex  breakdown  and  the  vortex  filament  method, 
offered  by  Cain  ( ref.  7),  have  given  much  insight  into  the  formative  stages  of  transient 
vortex  breakdown.  However,  in  attempt  to  better  model  the  flow  conditions  of 
axisymmetric  spiraling  flow  and  the  many  features  of  vortex  breakdown,  additional 
modifications  must  be  made  to  Cain’s  original  code.  This  section  will  present  these 
modifications  as  well  as  the  justification  and/or  rationale  for  their  employment. 

2.1  Upstream  Boundary  Conditions 

One  major  difficulty  with  any  numerical  simulation  lies  in  specifying  and  matching 
appropriate  boundary  conditions  within  the  numerical  domain.  These  difficulties  have 
been  well  documented  and  the  past  works  of  Nakamura  et  al  ( ref.  6),  Leonard  ( ref.  8),  and 
Garbowski  &  Berger  {ref.  9)  have  offered  guidance  for  many  subsequent  numerical 
methods. 

In  Cain’s  simulation,  all  flow  properties  along  the  boundary  of  the  flow  domain, 
with  the  exception  of  axial  vorticity  at  the  vorticity  source,  are  set  to  a  zero  value.  Such 
conditions  are  reasonable  for  the  radial  boundaries  conditions,  provided  that  the 
boundaries  are  a  significant  distance  away  from  vortex  core.  Computational  studies  by 
Brown  &  Lopez  {ref.  10)  and  experimental  studies  by  Harvey  {ref.  11)  concluded  that  a 
radial  boundary  should  be  placed  at  a  radius  of  at  least  two  core  diameters  from  the  center 
of  the  vortex  (centerline).  In  Cain’s  simulation,  this  boundary  is  approximately  four 
vortex  core  radii  from  the  centerline  and  therefore,  appropriate.  The  downstream 
boundary  condition,  similar  to  the  radial  boundary,  is  also  reasonable  so  long  as  the 
boundary  is  located  well  downstream  of  the  transient  vortex  filaments  and  their 
respective  induced  velocity  effect.  However,  the  upstream  boundary  condition,  in  this 
case  the  z  =  0  wall,  is  critical  in  calculating  the  velocity  and  vorticity  fields  within  the 
flow  domain.  Cain  elected  to  set  all  induced  flow  properties  at  the  wall  to  zero.  This 
condition  is  unsatisfactory  because  the  filaments  have  very  distinct  non- zero  velocities 
when  emerging  from  the  wall.  As  a  result,  the  effect  of  this  boundary  condition  taints  the 
interpolated  velocity  and  vorticity  fields  of  the  domain  near  the  upstream  wall  position. 


12 


This  inaccuracy  is  clearly  visible  in  Cain’s  vorticity  field  plots,  shown  in  Figure2-1  and 
Figure  2-2. 


0  0,5  1  1.5  2  25  3 
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Figure  2-1:  Azimuthal  vorticity  sign  switch,  meridional  plane,  T  =  1.2  ( ref.  7) 

Figure  2-1  illustrates  the  azimuthal  vorticity  field  of  the  domain.  Notice  that  all  along  the 
upstream  boundary,  z  =  0,  the  azimuthal  vorticity  is  zero.  As  previously  mentioned,  this 
inaccurate  condition  is  the  cause  of  the  discontinuity  in  the  flow  field  aft  of  the  upstream 
boundary,  z=0.  This  inaccuracy  is  also  illustrated  in  the  axial  vorticity  plot,  shown  in 
Figure  2-2.  At  the  time  shown  above  the  maximum  circulation  value  has  been  reach  and 
thus  the  vorticity  field  immediately  aft  of  the  wall  should  be  constant  with  time  (not 
radius)  and  clearly  non-zero. 
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Figure  2-2:  Axial  vorticity,  meridional  plane,  T  =  1.2  (ref.  7) 
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If  the  individual  properties  values  are  not  calculated  at  the  upstream  boundary 
condition,  as  in  Cain’s  case,  three  such  options  exist  to  correct  the  inaccurate  region, 
referred  to  as  Region  II.  The  first  would  require  placing  the  upstream  boundary  at  some 
location  downstream  of  the  z=0.  This  method  was  initially  utilized  by  Nakamura  et  al. 
{ref.  6)  in  their  past  numerical  simulations.  However,  Brown  &  Lopez  {ref.  10)  is  quick  to 
point  out  that  doing  such  requires  the  assumption  that  the  flow  at  the  upstream  boundary 
is  locally  independent  of  the  rest  of  the  flow  field  throughout  the  evolution  of  the  swirling 
flow.  Clearly  this  does  not  hold  true,  since  the  flow  field  in  this  region  varies  during  the 
circulation  and  velocity  ramp-up.  The  second  option  would  be  to  remove  or  ‘crop’  the 
inaccurate  portion  of  the  flow  field  from  the  data  range.  This  investigation,  however,  is 
interested  in  the  eventual  control  of  vortex  breakdown  within  the  flow.  These  controls 
may  include  the  ‘near-wall’  flow,  which  may  be  nestled  within  Region  II.  If  this 
investigation  were  concerned  with  only  the  far  downstream  structure  of  the  vortex,  this 
crop  option  would  be  suitable. 

The  third  option,  the  new  modified  numerical  simulation  allows  for  the 
calculations  of  the  wall  property  values  as  a  function  of  time;  specifically  those  of  the 
azimuthal  component  of  vorticity.  Therefore  the  inaccurate  region  can  be  eliminated  by 
applying  the  values  directly  to  wall  boundary  condition  at  each  time  step.  Velocity 
properties  at  the  wall  can  be  extrapolated  directly  from  the  data,  but  the  vorticity 
properties,  in  general,  and  the  azimuthal  vorticity  data  at  the  upstream  boundary  in 
particular,  requires  additional  calculations. 

The  Helmhotz  vortex  theorem  states  that  the  vorticity  moves  with  the  motion  of  a 
material  fluid  line.  This  asserts  that  as  the  induced  velocities  trigger  the  fluid  elements  to 
move,  they  additionally  move  the  vorticity  associated  with  the  fluid  elements.  This  does 
not  mean,  however,  that  the  vorticity  vector  and  velocity  vector  are  necessarily  aligned. 
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Figure  2-3:  Vortex  tube  orientation 

The  velocity  and  vorticity  vectors  can  be  decomposed  into  components  and  the  angles 
between  the  respective  vector  components  can  be  calculated.  These  angles,  €>u  and  ®0J, 
referred  to  as  the  helical  angles  of  the  velocity  and  vorticity,  respectively,  are  shown  in 
Figure  2-4. 


Figure  2-4:  Velocity  and  vorticity  vector  components 


To  be  consistent  with  past  studies,  this  investigation  will  utilize  the  a0  and  po  parameters, 
proposed  by  Brown  &  Lopez  ( ref.  10).  These  parameters  are  the  tangents  of  the  helix 
angle  for  the  velocity  and  vorticity,  respectively.  The  values  of  a  and  P  are  related  to  ®u 
and  3>w  by 
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a0  =  tan  (®u)  =  1 

uz„ 

Cl)  A 

Po  =  tan  (Ow)  =  — 


(2-1) 

(2-2) 


We  can  solve  for  the  tangent  of  the  helical  angle  of  velocity  at  the  wall,  a0,  using  the 
Equation  (2-3).  Note  that  a0  is  only  a  function  of  time  and  radial  position. 

ua  C  (t,r,z  =  0) 


a  =  tan  (Ou)  = 


u7 


U(t)  •  2  nr 


(2-3) 


where,  C  is  the  sum  of  all  circulation  contained  within  the  jth  concentric  ring  of  vortex 
filaments  and  U(t)  is  simply  the  free  stream  velocity.  As  previously  mentioned  the  total 
circulation  of  the  vortex  tube  is  discretized  into  nine  (9)  concentric  rings,  each  composed 
of  16  filaments,  using  a  Gaussian  distribution.  Therefore,  the  total  circulation  within  the 
jthring  is  simply  the  sum  of  all  the  distribution  constants,  G,  multiplied  by  the  total 
circulation,  Coo. 


C(t,rJ,z  =  0)  =  £2”'jc.r, 
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(2-4) 


Similarly,  since  we  are  modeling  the  case  of  unbounded  swirling  flow,  the  axial  vorticity 
at  z=0  is  defined  by  Equation  (2-5). 
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To  match  experimental  data,  we  again  choose  b  =  10.9  ( ref.  7).  Coo  is  specified,  as  it  is 
the  total  circulation  at  a  given  axial  location. 

Stream  surface  studies  by  Brown  &  Lopez  ( ref.  10)  have  shown  that  in  order  to 
produce  a  negative  azimuthal  vorticity  gradient,  the  ratio  of  a0  to  po  must  be  greater  than 
one.  Since  the  essential  feature  of  vortex  breakdown  is  the  generation  of  negative 

(X 

azimuthal  vorticity,  a  ratio  greater  than  one  must  be  used;  in  this  case  —of  1.5  is  chosen. 


Given  the  above  expression  for  the  axial  vorticity  at  the  boundary  condition,  coZ),  from 
(2-5),  a0  from  (2-3),  and  the  proportional  relationship  between  a0  and  po,  the  following 
expression  for  the  azimuthal  vorticity  at  the  wall  can  be  determined. 
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The  azimuthal  vorticity  is  a  function  of  both  time  and  radius.  Therefore,  the  upstream 
boundary  conditions  vary  depending  on  the  radial  distance  from  the  vortex  core,  as  well 
as  in  time.  Figure  (2-5)  shows  the  velocity  and  vorticity  distribution  along  the  modified 
vorticity  source  as  a  function  of  the  radial  distance  from  the  centerline.  Similarly, 

Figure  2-6  shows  the  “jet- like”  axial  velocity  profile  of  the  upstream  boundary  condition 
as  a  function  of  time. 


Figure  2-5:  Vortex  velocity  and  vorticity  distribution  at  z  =  0 

Figure  (2-5)  gives  insight  into  the  solid  body  and  irrotational  regions  of  the  vortex.  As 
with  all  real  vortices,  the  boundary  between  the  solid  body  region  and  irrotational  region 
occurs  at  the  radius  of  maximum  azimuthal  velocity.  This  radius  location  is  said  to  be  the 
edge  of  the  vortex  core.  Similar  to  Cain’s  simulation,  the  modified  simulation  has  a  core 
radius  size  of  0.30. 
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Axial  Velocity  as  a  Function  of  Radial  Distance 
Along  z  =  0  Plane 
x  =  0.2  sec  Cmax  =  20  Umax=5 


Figure  2-6:  Axial  velocity  profile  as  a  function  of  time  at  z  =  0 

Figure  2-6  illustrates  the  development  of  “jet-like”  axial  velocity  profile  at  the  upstream 
boundary  condition,  offered  by  the  modified  simulation.  At  T=0.9  the  profile  is 
established,  with  an  axial  velocity  at  the  center  of  the  vortex  core  of  almost  twice  the 
freestream  value.  It  is  important  to  note  that  steady  vortex  breakdown  LDA  experiments 
by  Escudier  ( ref.  12)  and  Faler  &  Leibovich  {ref.  13)  have  demonstrated  that  the  axial  jet 
profile  is  already  well  established  upstream  prior  to  breakdown.  However,  since  the 
simulation  models  transient  vortex  flow,  a  period  of  time  is  required  in  order  to  fully 
develop  this  profile. 

In  theory,  the  boundary  conditions  would  constitute  a  continuous  function. 
However,  in  order  to  implement  these  conditions,  numerically,  the  values  are  discretized 
radially  along  the  boundary  such  that  the  values  align  geometrically  with  the  outputted 
flow  domain  data.  The  benefit  of  such  alignment  is  that  no  additional  modifications  to 
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the  interpolation  subroutines,  structured  within  the  numerical  code,  are  needed.  An 
example  of  this  alignment  is  demonstrated  in  Figure  2-7. 


r 


Z=0  •  Outputted  nodal  data 

#  Discretized  boundary  data 


Figure  2-7:  Alignment  of  upstream  boundary  discretized  conditions 

The  hope  is  that  the  ability  to  better  capture  the  flow  immediately  aft  of  the  vorticity 
source  will  not  only  offer  further  insight  into  the  flow  mechanisms  but  will  present  a 
more  accurate  numerical  model. 


2.2  Convective  Velocity  Pitch-Up 

The  second  modification  to  Cain’s  original  code  is  the  introduction  of  a  non¬ 
constant,  time  dependent  convective  velocity  profile.  This  convective  velocity  is  defined 
as  the  velocity  at  which  the  fluid  mass  is  advected  through  the  domain  and  does  not 
include  induced  velocity.  Cain’s  original  code  fixed  the  axial  velocity,  at  which  the 
filaments  emerged  from  the  vorticity  source  and  moved  downstream,  to  a  constant  value. 
Simple  dimensional  analysis  can  show  that  circulation  is  simply  a  function  of  angle  of 
attack  (or  vane  angle),  velocity,  and  an  appropriate  length  scale,  whether  chord  length  or 
radius  of  the  vortex  core.  If,  in  transient  test,  the  angle  of  attack  and  length  scale  are  held 
constant,  such  as  in  delta  wing  or  bounded  pipe  flow  investigations,  circulation  then 
becomes  only  a  function  of  velocity.  With  this  in  mind,  the  goal  is  to  model  the 
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convective  velocity  proportionally  with  that  of  the  established  ‘pitch  up  and  hold’ 
circulation  profde. 

A  preliminary  attempt  to  alter  Cain’s  constant  ffee-stream  velocity  was  performed 
by  Wimer  {ref.  14).  Wimer  introduced  a  spatially  decreasing  convective  velocity  profile 
in  hopes  of  modeling  the  velocity  profile  associated  with  vortex  breakdown  over  a  delta 
wing.  The  axial  fluid  velocity  over  a  delta  wing  reaches  a  maximum  near  the  vorticity 
source  and  soon  then  rapidly  falls  towards  the  trailing  edge.  However,  Wimer  was  quick 
to  point  out  that  this  measured  axial  profile  over  such  a  delta  wing  is  a  result  of  steady 
state  condition,  a  condition  that  Cain  was  not  able  to  achieve.  Additionally,  he  notes  the 
difficultly  and  uncertainty  in  separating  the  convection  velocity  component  caused  by  the 
pressure  gradient  from  the  vortex  induced  velocity.  Despite  these  concerns,  he  enabled 
the  filaments  at  the  vorticity  source  to  convect  at  a  high  velocity  and  then  slow  down  with 
increasing  z-position. 

Unlike  Wimer’ s  spatially  dependent  convective  velocity  profile,  the  present 
modification  prescribes  the  convective  velocity  as  a  function  of  time,  which  enables  the 
convective  and  induced  velocities  to  remain  separate.  Under  this  condition,  all  flow 
within  the  domain  at  a  given  time  step  is  assigned  the  same  convective  velocity.  A 
qualitative  representation  of  this  condition  is  illustrated  in  Figure  (2-5).  This  condition 
maintains  continuity  and  better  models  the  flow  behavior  of  bounded  spiraling  flow,  such 
as  in  the  case  of  pipe  flow  {ref.  12). 


Z<0  Z=0  Z>0 

I 

'  Flow  Domain 


convective  [  convective  velocity  value 

velocity  profile  1  assigned  to  entire  flow  field 

Figure  2-8:  Convective  velocity  assignment  to  flow  field 
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Owing  to  the  aforementioned  proportional  relationship  between  velocity  and 
circulation,  a  baseline  velocity  profile  similar  to  the  circulation  profile  was  chosen,  which 
can  be  seen  in  Figure  (2-6).  Some  experimentation  was  done  in  arriving  at  the  present 
proportional  constant  of  one-fourth  (%).  In  keeping  with  Cain’s  color  designation 
scheme,  this  investigation  will  continue  to  use  the  color  red  to  indicate  a  region  of 
changing  circulation. 


Baseline  Simulation:  Circulation  and  Velocity  Profiles 


Figure  2-9:  Baseline  simulation:  circulation  and  velocity  profiles 

2.3  Induction  Turn-off  Zone 

Another  modification  made  to  Cain’s  original  simulation  is  the  creation  of  an 
induction  switch  off  zone.  For  reasons  which  will  be  further  explained  later  in  this 
investigation,  this  switch  allows  for  the  induction  on  the  vortex  filaments  within  a 
specified  region  , immediately  aft  of  the  vorticity  source,  to  be  turned  off.  This  turn-off 
means  that  the  entering  filaments  no  longer  receive  any  induced  velocities  from  the  other 
filaments  in  the  flow,  but  merely  convect  downstream  according  to  the  current  convective 
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velocity.  These  filaments,  however,  do  continue  to  induce  velocities  on  all  remaining 
filaments  downstream  of  the  turn-off  region.  As  the  simulation  progresses  and  new 
filaments  are  emitted  from  the  source,  they  too  are  turned  off  as  well.  Under  this 
condition,  the  turned-off  filaments  come  out  of  the  vorticity  source  axially  aligned  and 
continue  to  stay  straight  until  they  leave  the  induction  turn-off  region.  The  size  of  this 
turn-off  region  is  variable  and  user  specified.  This  switch  also  sets  the  stage  for  forced 
filament  bending.  In  this  case,  the  filaments  in  this  region  are  set  to  a  user-defined  tilt  or 
helical  angle.  This  scenario  will  be  later  employed. 

2.4  Removal  of  Blue  Region 

Cain’s  original  simulation  divides  the  vortex  tube  into  three  regions  -  green,  blue 
and  red.  The  green  shaded  filaments  indicate  regions  of  the  high,  final  circulation.  The 
blue  region  indicates  the  area  of  low  circulation  whereas  the  red  region  indicates  the  area 
that  circulation  is  increasing.  However,  preliminary  studies  by  Cain  (2001),  and  later 
studies  performed  by  Wimer,  showed  that  the  formative  features  of  breakdown,  including 
the  turning  point,  stagnation  point,  and  azimuthal  vorticity  sign  switch  location,  occur  in 
the  red  and  green  regions  only.  These  mechanisms  are  not  influenced  by  nor  affected  by 
the  nature  of  the  blue  region.  Therefore,  in  order  to  reduce  computational  run-time  and 
ease  memory  storage  requirements,  the  blue  region  of  constant  low  circulation  is 
removed.  As  a  result,  the  linear  velocity  and  circulation  pitch-up  begins  immediately  at 
the  start  of  the  simulation. 

2.5  Variable  Pitch-Up  Time 

One  of  the  key  features  studied  in  this  investigation  is  the  effect  of  the  velocity 
and  circulation  ramp-up  time  on  the  formative  features  and  locations  of  breakdown.  In 
order  to  do  such,  the  simulation  was  modified  to  allow  for  varying  ramp-up  time  values. 
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3  Simulation  Overview  &  Results 

Having  modified  the  vortex-filament  method  simulation,  this  section  will  present 
results  in  both  qualitative  and  quantitative  form.  Furthermore,  the  improvements  and 
differences  from  the  findings  offered  by  Cain  ( ref.  7)  will  be  addressed.  A  discussion  of 
the  simulation  results  is  also  presented  along  with  how  the  new  findings  provide  further 
evidence  for  the  inviscid  mechanism  of  the  self-induction  theory  of  vortex  breakdown. 
Finally,  this  section  will  also  offer  proof  that  the  modified  simulation  better  represents 
actual  vortex  flow  by  comparing  the  results  with  past  swirling  flow  experiments. 

3. 1  Vortex  Streakline  Representation 

The  solid  surface  rendition  will  be  presented  first  and  compared  with  Cain’s  results. 
In  these  results,  such  as  the  example  in  Figure  3-1,  the  vortex  lines  starting  initially  as  are 
rendered  as  surfaces  in  order  to  create  a  solid  shape  rather  than  a  tangled  assortment  of  all 
the  vortex  filaments,  seen  in  Figure  3-2.  These  vortex  lines  are,  in  fact,  streaklines  in 
that  they  represent  the  many  vortex  filaments  segments  released  successively  from  the 
fixed  vorticty  source,  at  one  instant  in  time.  Since  these  vortex  filaments  segments  are 
comparable  to  fluid  particles,  this  representation  offers  a  valid  comparison  to  dye  streaks 
used  in  past  flow  visualization  experiments. 


Figure  3-1:  Solid  representation,  r  =  0.08,  perspective  view 
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Streaklines,  Tau  :  0.20,  Time  :  1.75 


Figure  3-2:  All  vortex  filaments,  perspective  view 
The  only  slight  shortcoming  in  the  solid  representation  is  that  only  one  of  the  six  radii  can 
be  shown  at  one  particular  time.  This  restraint,  however,  is  comparable  to  flow 
visualization,  in  that  the  dye  gets  entrained  into,  and  remains  in,  a  single  radial  layer  of 
the  vortex  as  it  moves  with  the  fluid.  In  comparison,  this  minor  shortcoming  is  a  small 
price  to  pay  for  the  many  additional  intrinsic  benefits  in  visualizing  the  vortex  filaments. 

As  pointed  out  by  Cain  (ref  7),  one  such  benefit  to  the  vortex  filament  method  is 
the  ability  to  infer  induced  velocities  by  simply  looking  at  the  vortex  filaments 
orientation  relative  to  each  other.  Since  the  vorticity  vector  is  aligned  with  the  vortex 
filaments,  the  induced  velocity  from  any  one  particular  filament  on  the  flow  field  can  be 
determined  using  the  right  hand  rule.  This  principle  benefit  affords  insight  into  the 
formative  stages  of  vortex  breakdown  which  many  past  numerical  simulation  techniques 
could  not  offer. 

Three  different  sets  of  the  modified  simulation  results  are  plotted.  The  parameters 
of  the  baseline  simulation  are  shown  in  Table  3- 1 .  The  first  set,  shown  in  Figure  3-3, 
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shows  the  solid  representation  of  the  vortex  from  a  side  view.  The  vortex  filaments  in 
this  figure  emerge  from  the  vorticity  source  at  r  =  0.24. 

Table  3-1:  Baseline  simulation  parameters 


Parameter 

Value 

Units 

Initial  filament 

length,  Az 

0.1 

[length] 

Simulation  time 

step,  At 

0.01 

[time] 

Cinitial 

0.0 

[length]2  /  [time] 

C  final 

20 

[length]2  /  [time] 

Uinitial 

0.0 

[length]  /  [time] 

Ufinal 

5 

[length]  /  [time] 

^velocity 

0.20 

[time] 

^circulation 

0.20 

[time] 

Figure  3-3  also  shows  the  associated  radial  spokes.  These  radial  spokes  occur  in 
the  region  of  changing  circulation  and  are  needed  in  order  to  adhere  to  the  third  part  of 
the  Helmholtz  theorem,  which  requires  the  emergence  of  radial  vorticity  to  satisfy  the 
solenoidal  property.  Refer  to  Cain  {ref.  7)  for  a  more  detailed  explanation.  As  the  radial 
spokes  emerge  from  the  vorticity  source,  they  are  branched  off  with  circulation  equal  to 
AC  and  then  sent  all  the  way  far-field.  Although  the  branches  emerge  from  the  vorticity 
source  straight,  subsequently  they  do  not  remain  purely  in  the  radial  direction  due  to 
mutual  induction.  For  simplicity,  the  results  are  presented  in  time  intervals  of  t  =  0.25. 
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Figure  3-3:  Modified  simulation  results  with  radial  spokes,  r  =  0.24 
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Figure  3-3  (cont) 
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Streaklines,  Taj  :  0.20,  Time  :  2 


Figure  3-3  (cont) 

The  general  structure  and  evolution  of  the  vortex  in  our  modified  case  is  very  similar  to 
the  results  shown  by  Cain.  In  comparing  our  results  to  those  of  Cain,  it  is  important  to 
keep  in  mind  that  our  modified  simulation  does  not  include  the  low  circulation  blue 
region  used  by  Cain.  In  Cain’s  simulation,  the  ramp-up  region  (red  region)  emerged  from 
the  vorticity  source  immediately  after  the  blue  region,  at  T  =  0.6.  However,  with  no  blue 
region,  the  ramp-up  region  (red)  in  our  modified  simulation  begins  immediately  at  T  =  0. 
In  short,  the  start  of  our  modified  simulation  corresponds  to  Cain’s  simulation  at  T  =  0.6. 

The  second  set,  shown  in  Figure  3-4,  shows  the  same  view  of  the  solid 
representation  but  with  r  =  0.03,  the  inner  most  vortex  filaments. 


Streaklines,  Tau  :  0.20,  Time  :  0.25 
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Figure  3-4:  Modified  simulation  results,  r  =  0.03,  x  =  0.20 


27 


Streaklines,  Tau  :  0.20,  Time  :  0.75 


Streaklines,  Tau  :  0.20,  Time  :  1 


Streaklines,  Tau  :  0.20,  Time  :  1.25 


Figure  3-4  (cont.) 
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Streaklines,  Tau  :  0.20,  Time  :  2 
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Figure  3-4  (cont.) 

Unlike  the  r  =  0.24  results,  the  r  =  0.03  results  are  different  than  those  offered  by  Cain. 
The  early  development  of  the  vortex,  up  to  T  =  1.0,  matches  that  of  Cain’s  results  but 
thereafter  soon  differ.  In  Cain’s  simulation,  the  streaklines  of  the  inner  most  radii  do  not 
suddenly  explode  into  the  recirculation  bubble  seen  at  T=2.0;  instead,  the  behavior  of  the 
streaklines  is  less  dramatic  as  the  initial  tube  is  gradually  expanded  to  form  a  more 
streamline-shaped  recirculation  region.  A  side-by-side  comparison  of  this  is  illustrated 
below  in  Figure  3-5. 


M  odifi  ed  Simulation  Origi  rtal  Si  mulati  on  By  Cain 


Figure  3-5:  Comparison  of  inner  radius  streamlines  with  Cain 
The  third  set  of  plots  illustrates  the  interaction  between  the  inner  most  radius, 
r  =  0.03,  and  the  outer  most  radius,  r  =  0.42.  The  outside  filaments  are  rendered 
transparent  in  order  to  see  inside  of  the  vortex  tube  and  visualize  the  inner  most 
filaments,  shown  in  red.  It  should  be  noticed  that  the  both  the  inner  and  outer  radii 
filaments  initially  expanded  at  the  same  time  and  axial  location.  As  the  vortex  begin  to 
convect  downstream,  their  respective  expanding  regions,  or  bulges,  move  together. 
However,  as  the  vortex  tube  expands  further,  the  axial  velocity  of  the  inner  filaments 
decreases  due  to  the  increasing  presence  of  backflow.  This  decrease  in  velocity  causes 


29 


the  inner  filaments  to  lag  .  Furthermore,  as  the  backflow  increases,  the  inner  most 
filaments  wrap  around  the  opposing  reversed  flow  forming  the  core  of  recirculation 
bubble.  Since  it  is  primarily  the  initial  explosion  of  the  inner  core  which  feeds  the 
development  of  the  recirculation  bubble,  we  will  choose  to  view  the  inner  most  filaments 
in  determining  the  point  of  vortex  breakdown. 


Vortex  Filaments,  Tau  :  0.20,  Time  :  0.25 
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Figure  3-6:  Vortex  Filaments,  r  =  0.03,  r  =  0.42,  x  =  0.20 
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Vortex  Filaments,  Tau  :  0.20,  Time  :  2 


Figure  3-6  (cont.) 

The  behavior  of  the  streaklines  and  the  development  of  the  recirculation  bubble, 
at  the  inner  most  radius  in  particular,  as  offered  by  our  modified  simulation  models  the 
breakdown  phenomenon  better  than  Cain.  This  improvement  can  be  seen  when 
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comparing  our  streakline  results  to  that  of  past  flow  visualization  studies.  Such  a 
comparison  is  illustrated  in  Figure  3-7. 


Figure  3-7:  Vortex  simulation  comparison  to  experimental  investigation  {ref.  15) 
The  flow  visualization  pictures  in  the  top  half  of  Figure  3-7  are  from  a  transient  vortex 
breakdown  study  performed  by  Hamada  {ref.  15).  In  his  experiment,  a  needle  valve 
located  downstream  of  the  test  section  is  suddenly  opened  to  initiate  the  vortex 
breakdown  in  a  transient  manner.  It  is  this  transient  behavior  which  makes  the  results  by 
Hamada  an  ideal  source  for  comparison.  The  simulation  results  match  well  with  the 
flow  visualization.  One  slight  difference,  though,  is  that  the  outer  edge  of  both  the  initial 
and  final  bubble  in  the  flow  visualization  is  smoother  than  those  of  the  simulation.  It  is 
important  to  remember  that  this  is  an  inviscid  simulation,  and  the  smoothing  effects  of 
viscosity  are  not  present.  This  inviscid  nature  is  the  reason  for  the  jaggedness  of  the 
simulated  bubbles.  Similarly,  it  is  the  lack  of  viscosity  which  does  not  allow  the 
recirculation  bubble  to  reach  a  final  stationary  position  in  our  simulation.  In  real  flow, 
filaments  with  equal  but  opposite  vorticity  would  cancel  due  to  viscosity  effects. 
However  the  simulation  is  not  able  to  achieve  this  annihilating  effect,  and  as  a  result, 
continues  to  impose  undiminished  induction  from  the  filaments  on  the  rest  of  the  field; 


32 


this  continuing  mutual  induction  of  filaments  from  within  the  recirculation  bubble 
causes  the  bubble  to  convect  downstream.  Conversely,  in  real  flow,  the  smoothing  effect 
of  viscosity  would  cancel  out  inflowing  filaments,  allowing  the  bubble  to  attain 
equilibrium  position  in  the  final  state. 

3.2  Velocity  and  Vorticity  Fields 

In  this  section,  we  present  the  azimuthal  vorticity  and  radial  and  axial  velocity  as 
calculated  along  the  meridional  plane  of  the  vortex  tube.  To  aid  visualizing  the 
instantaneous  vortex  shape,  vortex  filaments  at  r  =  0.12  and  r  =  0.42  are  superposed  onto 
the  plots.  All  three  sets  of  plots  show  much  improvement  over  Cain’s  results 
immediately  aft  of  the  vorticity  source,  z  =  0  boundary.  One  example  of  this 
improvement  can  be  in  the  radial  velocity  plots.  Cain’s  results  showed  a  relatively  large 
region  of  radial  contraction  immediately  aft  of  the  vorticity  source.  We  will  refer  to  this 
contraction  behavior  as  the  “necking  effect.”  The  modified  simulation,  with  its  improved 
upstream  boundary  conditions,  is  able  to  greatly  reduce  the  necking  effect,  both  in 
magnitude  and  in  duration.  Another  improvement  can  be  seen  in  the  plots  of  azimuthal 
vorticity.  Similar  to  the  radial  velocity  case,  the  modified  simulation  is  able  to  correct  for 
inaccuracy  of  Cain’s  results  immediately  aft  of  the  vorticity  source;  the  azimuthal 
vorticity  no  longer  approaches  a  zero  value  at  the  upstream  boundary.  These 
improvements  are  not  only  required  to  better  model  vortex  behavior  but  are  critical  to  the 
projected  streamline  study,  which  will  be  presented  in  the  subsequent  chapters. 

3.2. 1  Color  contour  azimuthal  vorticity  plots 

Azimuthal  Vorticity,  Meridional  Plane,  Time  =  0.4 
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Figure  3-8:  Azimuthal  vorticity  color  contour  plots 
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Azimuthal  Vorticity,  Meridional  Plane,  Time  =  1.6 


Figure  3-8  (cont.) 
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Azimuthal  Vorticity,  Meridional  Plane,  Time  =  2 
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Figure  3-8  (cont.) 


3.2.2  Color  contour  radial  velocity  plots 

1.5 
1 

0.5 

-  0 
-0.5 
-1 
-1.5 

1  23456789  10 


Radial  Velocity,  Meridional  Plane,  Time 


Figure  3-9:  Radial  velocity  color  contour  plots 
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Radial  Velocity,  Meridional  Plane,  Time  =  2 
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Figure  3-9  (cont.) 
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3.2. 3  Color  contour  axial  velocity  plots 


Axial  Velocity,  Meridional  Plane,  Time  =  0.4 
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Axial  Velocity,  Meridional  Plane,  Time  =  1.2 
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Figure  3-10:  Axial  velocity  color  contour  plots 
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Axial  Velocity,  Meridional  Plane,  Time  =  2 


Figure  3-10  (cont.) 


3.3  Validation  of  Breakdown  Features 

Past  investigations  have  highlighted  four  principle  flow  features  associated  with 
breakdown:  axial  pile-up,  radial  expansion  and  contraction,  a  vorticity  sign-switch  and  a 
flow  turning  point.  Cain’s  simulation  gave  evidence  to  the  argument  that  all  four  of  these 
features  are  interlinked  result  of  the  self-induced,  inviscid  feedback  mechanism  of  vortex 
breakdown.  This  section  offers  brief  explanations  of  the  four  properties  and  gives 
evidence  to  their  respective  manifestation  in  the  modified  simulation  results. 

3.3.1  Radial  expansion 

Past  studies  have  clearly  shown  that  radial  expansion  arises  directly  from  an 
azimuthal  vorticity  gradient.  In  order  to  offer  an  explanation  for  radial  expansion,  it  is 
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best  to  present  this  azimuthal  vorticity  as  a  meridional  cut  of  discretized  vortices.  This 
meridional  cut  is  illustrated  in  Figure  3-11.  The  azimuthal  vorticity  of  equal  strength 
upstream  decreases  towards  downstream  and  no  vorticity  exists  beyond.  At  the  first 
vortex,  where  there  is  no  preceding  vortex  in  the  downstream  direction,  all  the  induction 
are  additive  such  that  the  resultant  is  a  radially  outward  velocity.  As  you  move  upstream, 
where  the  vorticity  is  of  equal  magnitude,  the  radially  outward  and  inward  induction  by 
neighboring  vortices  cancel  out.  Extending  this  to  a  3-D  case,  the  induced  radial 
velocity,  produced  by  the  azimuthal  vorticity  gradient,  causes  the  vortex  tube  to  expand 
in  the  downstream  direction;  once  the  azimuthal  disappears,  the  radial  size  of  the  vortex 
remains  the  same,  as  seen  in  the  radial  velocity  plots  in  chapter  3.2.3  near  z=0. 


decreasing  azimuthal  coq 


downstream 


upstream 


Figure  3-11:  Radial  expansion  ( ref.  7) 


3.3.2  Axial  pile-up 

The  effect  of  the  radial  expansion,  discussed  in  the  previous  section,  also  leads  to 
axial  pile-up.  As  noted  by  Cain,  pile-up  refers  to  the  axial  contraction  of  the  vortex  tube. 
Similar  to  radial  expansion,  pile-up  is  also  a  result  of  an  azimuthal  vorticity  gradient.  As 
with  the  previous  section,  we  can  use  a  meridional  cut  of  discretized  vortices,  shown  in 
Figure  3-12,  to  illustrate  the  creation  of  pile-up. 
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Figure  3-12:  Axial  pile-up  (ref.  7) 

Since  the  axially  induced  velocity  at  the  center  of  the  vortex  is  induced  symmetrically, 
we  can  further  simplify  our  example  by  considering  the  effects  of  only  the  top  or  bottom 
half  row  of  vortices.  In  this  case,  we  will  look  at  the  top  row.  The  induced  axial  velocity 
at  point  (a)  is  generated  only  from  all  the  upstream  vortices;  however,  the  greater  part  of 
the  contribution  comes  from  only  a  few  neighboring  upstream  vortices.  As  you  move 
further  upstream,  at  point  (c),  the  contributing  pool  of  neighboring  vortices  include  those 
downstream  as  well  as  upstream.  As  a  result,  induction  comes  from  not  only  vortex  (3) 
but  from  (2)  and  (1)  as  well.  The  additive  induced  velocities  from  both  (2)  and  (1)  cause 
the  net  induced  axial  velocity  at  (c)  to  be  greater  than  either  (b)  and  (a).  As  you  move 
upstream  and  further  away  from  the  effects  of  vortex  (1)  and  vortex  (2),  the  induced  axial 
velocity  approaches  a  constant  value.  In  the  case  of  our  simulation,  this  constant 
condition  is  established  as  the  region  of  circulation  change  moves  well  downstream. 

Such  constant  axial  velocity  is  illustrated  in  the  previous  contour  plots  of  axial  velocity. 
Observe  also  that  the  radial  profile  of  the  axial  velocity  reaches  maximum  at  centerline,  a 
jet-like  profile.  This  is  also  analogous  to  the  jet  profile  which  develops  at  the  leading 
edge  of  a  delta  wing  vortex.  Contrary  to  this  case,  a  positive  vorticity  gradient  would 
cause  the  induced  velocities  to  be  reversed.  These  reversed  radial  velocities  would  lead 
to  axial  acceleration  and  a  contraction  at  the  downstream  edge. 

The  axial  velocity  gradient,  illustrated  in  Figure  3-12,  causes  the  axial 
components  of  the  vortex  filaments  to  contract.  This  axial  contraction  leads  to  the 
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increase  in  the  cross  sectional  of  the  vortex  segment.  Since  the  local  circulation 
around  the  vortex  segment  remains  constant,  this  leads  to  decrease  vorticity,  thereafter,  to 
a  corresponding  decrease  in  azimuthal  vorticity.  This  decrease  in  turns  creates  a  stronger 
azimuthal  vorticity  gradient  and  more  pile-up.  Therefore,  once  pile-up  begins,  this 
endless  feedback  cycle  causes  the  effect  of  pile-up  to  dramatically  increase  with  time. 

3.3.3  Vorticity  sign-switch 

As  previously  noted  in  the  explanation  of  axial  pile-up,  there  is  a  decrease  in 
vorticity  in  the  location  of  a  vorticity  gradient.  Hence,  we  would  similarly  expect  a 
decrease  in  azimuthal  vorticity  in  the  regions  of  radial  expansion  and  pile-up.  Cain 
(2001)  also  noted  that  the  azimuthal  velocity  gradient  causes  tilting  in  the  vortex 
filaments  in  the  region  of  radial  expansion  and  pile-up.  Due  to  this  tilting  behavior,  the 
vector  decomposition  of  the  vortex  filament  will  illustrate  a  reverse  in  azimuthal 
vorticity,  referred  to  as  the  sign-switch.  This  tilting  behavior  is  illustrated  in  Figure  3-13. 


Figure  3-13:  Azimuthal  vorticity  sign-switch  ( ref.  7) 

As  you  can  see  in  Figure  3-13,  the  azimuthal  velocity  at  the  downstream  location  has 
decreased  from  that  of  the  upstream  location  due  to  the  radial  expansion  and  subsequent 
axial  pile-up.  If  we  then  decompose  the  vector  components  of  the  vortex  filament,  we 
soon  uncover  the  azimuthal  vorticity  sign-switch.  Again,  we  can  decompose  the  vorticity 


vector  because  of  the  fact  that  the  vector  is  always  aligned  with  the  vortex  filament. 
This  vortex  tilting  is  furthered  observed  in  our  simulation  results,  as  shown  in 
Figure  3-14. 
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Vortex  Filaments,  Time:  0.5 


Figure  3-14:  Azimuthal  vorticity  sign-switch,  simulation  results 
Figure  3-14  illustrates  the  vortex  filaments  at  T  =  0.5  as  well  as  the  vorticity  components 
projected  along  the  meridional  plane. 

3.3.4  Turning  point 

The  turning  point  is  the  location  where  the  vortex  filaments  have  turned  inward 
on  themselves.  This  culminating  location  is  marked  by  a  region  of  radial  contraction  and 
reversed  axial  flow  inside  of  the  vortex  tube  (backflow).  These  two  properties  are  both 
resultants  of  the  induced  velocities  triggered  by  the  azimuthal  sign  switch.  The  turning 
point  can  further  be  defined  as  the  location  where  the  axial  vorticity  is  zero,  such  that  the 
axial  vorticity  upstream  of  this  point  is  positive  and  the  velocity  downstream  along  the 
filament  is  negative.  This  distinct  point  of  separation  manifests  itself,  in  the  vortex 
filament  simulation,  as  the  point  where  the  vortex  filament  is  vertical  relative  to  the 


z-axis. 
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3.3.5  Property  snap-shot 

As  supported  by  Cain  (2001),  vortex  breakdown  can  be  initiated  and  defined 
through  the  entirely  inviscid  mechanism  of  pile-up,  radial  expansion  and  contraction,  and 
vorticity  sign  switching.  The  intertwining  of  these  three  behaviors,  all  of  which  are 
caused  by  a  simple  azimuthal  vorticty  gradient,  culminates  in  the  turning  point.  Once  the 
turning  point  has  been  formed,  the  downstream  region  forms  a  dissipative,  turbulent  flow 
that  leads  to  a  steady  state  vortex  breakdown.  Similarly,  the  flow  immediately  upstream 
of  the  turning  point  continues  to  turn  in  on  itself,  creating  a  recirculation  region.  Within 
this  region,  the  counter  rotating  flow  further  interacts  with  the  upstream  flow,  leading  to 
viscous  activity.  It  is  the  dissipative  viscous  activity,  along  with  the  reversed  axial  flow, 
which  leads  to  both  the  stagnation  of  the  recirculation  region  and  establishment  of  steady 
vortex  breakdown.  Visualizations  of  the  dissipative  viscous  activity  within  a  nearly 


axisymmetric  recirculation  region  are  offered  by  Sarpkaya  ( ref.  16)  in  the  Figure  3-15. 


Figure  3-15:  Viscous  activity  within  an  unsteady  recirculation  region  (ref.  15) 

In  summary,  the  results  of  the  vortex  fdament  simulation  support  the  notion  that  transient 
formation  of  vortex  breakdown  is  not  only  an  inviscid  self-induced  process  but  initiated 
by  an  azimuthal  vorticity  gradient. 
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The  vortex  filament  method  is  a  very  powerful  method  as  it  enables  us  to  see 
all  four  features  of  the  inviscid  self-induction  mechanism  in  one  simple  snapshot  of  time. 
It  is  especially  useful  in  illuminating  the  movement  of  the  vorticity  vector  which  soon 
manifests  itself  as  the  sign  switch.  One  such  snapshot  is  illustrated  in  Figure  3-16. 


Figure  3-16:  Properties  of  the  inviscid  self-induction  mechanism 


3. 4  Swirl  Number  In  vestigation 

Many  non-dimensional  parameters  are  considered  and  associated  with  the 
characterization  of  vortex  flows  have  been  considered.  Dimensional  considerations  show 
that  a  global  description  of  vortex  flow  must  include  some  equivalent  swirl  parameter. 

As  pointed  out  by  Escudier  {ref.  12),  a  swirl  parameter  can  be  defined  in  terms  of 
numerous  possible  flow  properties,  to  include  axial  flux  and  angular  momentum  flux. 
However,  one  such  common  and  simple  choice  to  define  a  swirl  parameter  is  in  terms  of 
the  bulk  vortex  velocity.  By  means  of  this  velocity,  the  swirl  parameter  is  expressed  as 


wherein  Coo  is  the  total  circulation,  Uoo  is  the  free  stream  velocity  and  L  represents  some 
length  scale.  This  particular  swirl  parameter  is  commonly  referred  to  as  the  swirl 
number,  and  has  been  used  in  many  guidevane  induced  swirl  studies.  Since  Coo,  which 
has  the  dimension  of  velocity  and  length,  must  depend  on  the  angle  of  attack,  or  guide 
van  angle,  a ,  Q  is  uniquely  determined  by  a. 

With  this  swirl  number,  this  investigation  is  able  to  make  a  connection  between, 
and  further  compare,  the  results  of  the  vortex-filament  method  to  those  of  experimental 
findings  offered  by  Wimer  {ref.  14).  In  his  digital  particle  image  velocimetry  (DPIV) 
investigation  of  vortex  breakdown,  Wimer  quantified  flow  fields  of  breakdown  over  a 
delta  wing.  With  his  measurements,  he  was  able  to  calculate  the  circulation  and  velocity 
fields  of  the  vortex  flow  over  a  delta  wing  at  various  angles  of  attack  (AO A),  in  relation 
to  the  location  of  breakdown.  By  specifying  the  swirl  number,  Q,  in  terms  of  the  present 
simulation  input  parameters,  and  equating  this  Q  to  corresponding  experimental  value  we 
can  associate  the  results  of  the  vortex-filament  method  with  his  data  at  given  a  angle  of 
attack.  Figure  3-17  illustrates  vortex  streaklines,  of  r  =  0.03,  for  various  angles  of  attack. 
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Figure  3-17:  Simulation  results  for  various  equivalent  angles  of  attack,  r=0.03 
As  you  can  see  from  Figure  3-17,  vortex  breakdown  has  only  occurred  at  the  30  deg 
angle  of  attack  at  T=1.7.  The  3  deg  vortex  flow  has  translated  furthest  downstream 
because  it  does  not  experience  the  rigorous  reverse  flow  phenomenon  which  is  evident  in 
the  other  two  flows.  Although  the  vortex  flow  in  the  15  deg  and  3  deg  angle  of  attack 
case  shows  signs  of  radial  expansion  and  possible  reversed  flow,  no  turning  point  has 
been  observed.  In  fact,  as  time  increases,  the  flow  in  both  of  these  two  cases  continues  to 
expand  and  transient  downstream,  but  never  reaches  a  stage  to  form  a  turning  point. 
Although,  we  can  not  argue  that  the  flow  never  breaks  down,  but  rather  can  argue  that  the 
flow  never  breaks  down  within  our  established  simulated  domain.  If  we  think  of  our 
domain  in  terms  of  a  wing  on  a  aircraft,  the  behavior  represented  in  the  3  deg  and  15  deg 
case  would  be  analogous  to  vortex  breakdown  that  occurs  aft  of  the  trailing  edge.  Our 
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findings  hold  true  with  Wimer,  who  noted  that  at  angles  less  than  15  deg  illustrated  no 
breakdown  occur  over  the  delta  wing  surface 

A  point  to  note  is  that  the  delta  wing  used  by  Wimer  has  a  finite  chord  length, 
whereas  the  only  length  scale  in  our  simulation  is  that  of  the  vortex  core  radius. 

However,  we  can  use  the  30  deg  case  as  a  baseline  case  to  define  the  effective  simulation 
chord  length,  ceff,  and  by  equating  the  swirl  number  of  Wimer 

q  _  C  gym  ax 

Wimer  Uqo  c 

to  that  of  the  present 

c 

o  _  oofinal 

“VFM  "  tj  c 

u°o  ceff 

and  this  ceff,  is  used  for  other  angles  of  attack.  With  this  method  we  can  further  compare 
the  breakdown  locations  of  other  equivalent  simulations  to  those  recorded  by  Wimer 
(i ref.  14)  and  Srigrarom  ( ref.  17).  The  simulation  for  various  swirl  numbers  is  carried  out 
by  only  changing  the  Cfmai  value  of  Table  3-1,  with  the  other  parameters  remaining 
unchanged.  For  means  of  comparison,  the  simulated  location  of  breakdown,  Xb  ,  is 
defined  to  be  the  location  at  which  the  recirculation  bubble  has  reached  95%  of  its 
maximum  radial  size.  These  comparisons  are  shown  in  Table  3-2. 


Table  3-2:  Swirl  Number  (Q)  vs.  experimental  and  simulated  breakdown  locations 


Swirl 

Number 

Equivalent 

Angle  of 

Attack 

Simulated 

Breakdown 

Location,  Xb 

Xb  /  Ceff 

(VFM) 

Xb/c 

(Wimer) 

Xb/c 

(Srigrarom) 

Q=  13.33 

~30  deg 

3.5 

30% 

30% 

-30% 

Q=  11.11 

~25  deg 

5.2 

45% 

48% 

-40% 

Q.  =  8.88 

~20  deg 

6.5 

56% 

60% 

-70% 

Q  =  6.66 

~15  deg 

None 

- 

None 

None 

Q=  1.33 

~3  deg 

None 

- 

None 

None 
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Since  our  simulation  models  an  unsteady  process,  while  steady-state  is  reached  in 
experiments,  the  above  method  gives  only  a  rough  approximation  of  simulated  vortex 
breakdown  location.  However,  regardless  of  this  fact,  the  data  offered  in  Table  3-2 
points  out  the  effectiveness  of  the  simulation  to  determine  breakdown  locations  in 
reference  to  the  experimentally-matched  30  deg  baseline  case. 

One  noteworthy  criterion  is  the  minimum  swirl  number  value  at  which  vortex  breakdown 
first  occurs  within  our  effective  domain.  This  threshold,  referred  to  as  the  critical  swirl 
number,  ^critical ,  was  obtained  by  slowly  narrowing  down  the  swirl  number  range 
between  Q  =  8.88  and  Q  =  6.66.  A  critical  swirl  number  of  approximately  7.1  was 
reached.  This  critical  value  is  equivalent  to  a  65  deg  swept  delta  wing  set  at  a  1 7  deg 
angle  of  attack. 

3.5  Ramp-up  (time)  Investigation 

One  of  the  added  benefits  of  the  modified  vortex- filament  method  is  the  ability  to 
investigate  the  formative  stages  of  breakdown  at  varying  ramp-up  times.  The  ramp-up 
time  is  the  time  required  for  the  circulation  and  velocity  profile  to  change  from  a 
condition  at  rest  (C=0,  U=0)  to  the  maximum  user-defined  values.  In  physical  terms,  this 
ramp-up  time  is  analogous  to  the  time  required  for  a  wing  to  pitch-up  from  level  flight  to 
a  specified  angle  of  attack.  All  results  up  to  this  point  have  used  a  ramp-up  time  of 
x  =  0.20.  From  experiments,  we  know  that  the  resultant  shape  and  orientation  of  bubble 
breakdown  should  be  independent  of  pitch-up  time,  so  long  as  the  same  final  swirl 
number  value  is  reached.  In  physical  terms,  the  breakdown  over  a  wing  does  not  matter 
how  long  it  takes  for  the  wing  to  pitch  up  to  a  specified  angle  of  attack,  so  long  as  the 
same  final  angle  of  attack  is  reached.  Therefore,  we  would  expect  no  differences  in  the 
formative  stages  of  breakdown  with  varying  ramp-up  times;  additionally,  we  would  also 
expect  to  see  all  four  previously  discussed  features  of  the  inviscid  self-induction 
mechanism.  Figure  3-18  shows  the  modified  simulation  results  of  r=0.24  and  r  =  0.03  for 
the  Q  =  13.33  case  at  various  ramp-up  times. 

One  notable  feature  of  Figure  3-18,  is  the  difference  in  breakdown  locations.  The 
quicker  the  ramp-up  (smaller  the  value  of  x),  the  faster  breakdown  occurs.  Likewise,  the 
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quicker  the  ramp-up,  the  quicker  the  initial  onset  of  radial  expansion  and  azimuthal 
vorticity  sign-switch. 


Vortex  Streaklines,  R  =  0.24,  0.03,  Time:  1 .03,  Tau:  0.005 
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Figure  3-18:  Vortex  streaklines  at  breakdown  for  varying  ramp-up  times 
In  general,  the  smaller  the  x  value,  the  faster  all  aforementioned  features  of  the  inviscid 
self-induction  mechanism  are  observed.  This  trend  is  attributed  to  the  differential  rate  of 
change  in  circulation.  A  faster  differential  rate  of  change  in  circulation  causes  a  stronger 
vorticity  gradient,  which  ultimately  leads  to  faster  breakdown.  This  trend  is  comparable 
to  that  of  wingtip  vortices  on  an  aircraft.  As  an  aircraft  pitches  up,  it  develops  more  lift 
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over  the  wings.  This  lift  force  is  directly  related  to  the  circulation  around  the  wing 
section.  Therefore,  the  faster  the  pitch  up,  the  larger  the  differential  rate  of  change  in 
circulation  over  the  wing  and  the  faster/tighter  the  roll-up  of  the  wingtip  vortices.  In  an 
effort  to  relate  the  onset  of  radial  expansion  to  the  formation  of  the  breakdown  bubble, 
we  will  introduce  a  new  non-dimensionalized  parameter,  T*.  In  this  case,  T*  is  the  ratio 
of  the  time  at  which  the  breakdown  bubble  is  95%  expanded,  Tb,  to  the  initial  radial 
expansion  time  of  the  vortex  core,  Te. 


y*  _  Tb 


The  values  of  TB,  TE  and  T*  for  varying  ramp-up  times  are  shown  below  in  Table  3-3. 
Please  note  that  all  three  cases  represented  in  the  table  were  run  with  Cfmai  =  20. 


Table  3-3:  Non-dimensional  time  (T*)  vs.  ramp-up  time  (x) 


Ramp-up 

time,  x 

Te 

Tb 

x  =  0.20 

0.21 

1.04 

4.95 

x  =  0.05 

0.28 

1.39 

4.97 

x  =  0.005 

0.35 

1.72 

4.91 

As  you  can  see  from  the  table,  although  the  values  of  TE  and  TB  differ  for  each  case,  the 
values  of  T*  are  very  similar.  This  fact  would  lead  us  to  argue  that  there  is  a  unique 
relationship  between  the  time  of  initial  expansion  and  the  time  of  breakdown. 
Furthermore,  it  suggests  that  once  the  inviscid  self-induction  mechanism  causes  radial 
expansion,  the  subsequent  process  of  breakdown  is  already  predetermined.  If  this  is  true, 
then  any  attempt  to  control  or  manipulate  breakdown  must  be  focused  not  only  at  the  far 
upstream  position  but  early  in  the  formative  stages  of  breakdown  -  prior  to  the  time  of 
initial  expansion.  However,  in  order  to  further  validate  this  argument,  further 
simulations,  at  varying  Cfinai  values,  and  experiments  need  to  be  investigated. 
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4  Qualitative  Investigation 

The  focus  of  this  section  will  be  to  present  a  qualitative  study  of  the  shock-like 
behaviors  exhibited  in  the  preceding  vortex  filament  simulation  results.  A  discussion  of 
these  results  is  offered  and  compared  with  the  traditional  theory  that  breakdown  behaves 
much  like  the  sudden  transition  observed  in  a  shock  wave  or  hydraulic  jump. 

Additionally,  this  section  also  gives  insight  into  the  scenario  for  the  generation  and 
subsequent  growth  of  vorticity  gradients,  as  well  as  their  influence  on  radial  expansion. 

4.1  Elements  of  Shock-like  Behavior 

Numerous  theoretical  investigations  with  the  aim  of  grasping  the  complex 
phenomenon  of  vortex  breakdown  have  been  made  since  the  first  reported  observation  of 
breakdown.  Among  the  vast  literature,  many  studies,  Benjamin  ( ref.  18)  and  Keller  et  al 
( ref.  19)  in  particular,  have  considered  the  vortex  breakdown  problem  from  a  shock-like 
perspective,  making  the  analogy  between  a  hydraulic  jump  and  steady  vortex  breakdown. 
Although  the  vortex-filament  method  does  not  model  the  wave  effects  as  such,  the 
methods  does  provide  valuable  insight  into  the  temporal  evolution  of  flow  parameters. 

For  example,  the  gradual  build-up  of  the  axial  velocity  and  vorticity,  both  of  which 
commence  at  the  onset  of  the  simulation,  culminate  into  a  sudden  jump  or  switch  in 
magnitude.  These  sudden  jumps  in  flow  properties  are  manifested  as  the  axial  velocity 
stagnation  location  and  sudden  vorticity  sign  switch;  both  are  present  during  breakdown. 

4.1.1  Axial  Pile-up 

In  addition  to  the  raw  velocity  and  vorticity  data,  the  vortex  filament  method  also 
enables  us  to  look  at  the  transient  evolution  of  the  filaments,  subject  to  induced  velocities, 
in  their  relation  to  each  other.  Since  the  filaments  can  become  stretched  and/or  shortened 
during  the  simulation,  it  is  best  to  visualize  the  axial  location  of  the  filament  nodes  rather 
than  the  actual  filaments  themselves.  Again,  it  is  important  to  keep  in  mind  that  our 
simulation  models  the  formative  stages  of  transient  breakdown.  Figure  4-1  illustrates  the 
axial  position  of  nodes,  pertaining  to  the  inner-most  radius,  as  a  function  of  time.  To  help 
in  visualizing  the  relative  movement  of  the  nodes,  dotted  lines  are  added  connecting 
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identical  nodes.  In  other  words,  the  dotted  lines  represent  axial  trajectories  of  specific 
nodes.  As  you  can  see,  the  nodes  do  not  maintain  the  initial  equi-distant  spacing  (dz  = 
0.1)  as  the  filaments  begin  to  lengthen  and  shorten.  This  variation  in  nodal  spacing  is 
attributed  to  the  difference  in  induced  velocities  caused  by  the  self-induction  process. 


Axial  Node  Pile-Up  Investigation 
Similar  Linear  Velocity  and  Circulation  Ramp-Up 


Figure  4-1 :  Shock  behavior:  axial  node  pile-up,  short  range 

As  time  progresses,  the  spacing  of  incoming  nodes  is  lengthened  while  the  spacing 
between  the  nodes  further  downstream  begins  to  shorten.  This  compacting  behavior  is 
referred  to  as  ‘pile-up’  and  is  a  result  of  the  decreasing  axial  velocity  downstream.  This 
decreasing  axial  velocity  is  a  precursor  to  the  reverse  flow  seen  immediately  prior  to 
vortex  breakdown. 

While  Figure  4-1  illustrates  the  transitive  effect  of  pile-up  just  prior  to  vortex 
breakdown,  Figure  4-2  illustrates  the  continuing  transitive  effect  after  vortex  breakdown 
has  occurred.  In  this  case,  the  number  of  nodes  within  the  pile-up  region  has  greatly 
increased,  even  though  the  axial  width  of  the  region  has  not.  Clearly,  this  pile-up  effect 
is  analogous  to  the  temporal  formation  of  a  shockwave.  Since  our  simulation  is  a 
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transient  study,  the  recirculation  region  and  stagnation  point  is  continuously  convected 
downstream.  It  is  because  of  this  reason  that  the  pile-up  region,  which  is  closely  located 
to  the  recirculation  region,  continues  to  transient  downstream  in  Figure  4-2.  Similar  to 
Figure  4-1,  lines  are  used  to  trace  the  relative  movement  of  similar  nodes.  However,  in 
this  case,  the  colored  lines  are  shown  to  highlight  an  important  behavior  -  nodal 
overlapping.  If  you  look  at  T=1 .5,  you  will  see  the  green,  blue  and  red  lines,  in  ascending 
order.  However,  at  T=2.0,  the  order  of  the  lines  are  reversed.  This  inversion  is  caused  as 
downstream  nodes  are  pulled  or  “sucked”  into  the  recirculation  region,  overlapping  the 
nodes  immediately  behind  them.  This  continuing  overlapping  behavior  can  also  been 
seen  experimentally,  as  flow  passing  over  the  breakdown  shear  layer  is  reversed  and 
circulated  back  through  the  bubble.  It  is  important  to  note  that  the  simulation  results  give 
rise  to  a  systematic  overlapping  behavior  during  breakdown.  By  systematic,  we  mean 
that  the  behavior  follows  an  orderly  sequence  and  is  not  simply  a  random  scattering  of 
nodes. 

Axial  Node  Pile-Up  Investigation 
Similar  Linear  Velocity  and  Circulation  Ramp-Up 


Figure  4-2:  Shock  behavior:  axial  node  pile-up,  long  range 
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Another  method  of  illustrating  the  shock-like  behavior  of  vortex  breakdown 
includes  graphing  the  inner-radius  (r=0.03)  nodal  spacing  between  adjacent  nodes,  Az,  as 
a  function  of  axial  distance.  This  method  is  shown  in  Figure  4-3  and  Figure  4-4.  Similar 
to  the  previous  plots,  Figure  4-3  also  illustrates  a  transient  shock-like  behavior.  By 
T  =  0.8,  which  is  just  prior  to  breakdown,  the  spacing  of  the  vast  majority  of  nodes  is 
much  less  than  the  initial  dz=0. 1  value.  In  fact,  the  downstream  nodal  spacing  appears  to 
approach  Az=0  as  time  progresses.  The  negative  Az  values  represent  the  nodal 
overlapping  effect.  According  to  the  figure,  nodal  overlapping  begins  at  approximately 
T=1.6.  This  time  matches  well  with  the  previous  plots,  as  it  is  exactly  the  same  time  we 
saw  the  colored  lines  begin  to  overlap,  in  Figure  4-2. 


Axial  Node  Pile-Up  Investigation 
Similar  Linear  Velocity  and  Circulation  Ramp-Up 
x  =  0.2  sec  Cmax  =  20  Umax  =  5 


Figure  4-3:  Nodal  spacing  (Az)  vs.  axial  distance,  long  range 

Although  Figure  4-3  shows  a  long  range  time  history  of  the  transient  shock-like  behavior 
and  axial  pile-up,  it  does  a  poor  job  at  showing  the  early  transient  behavior.  In  contrast, 
Figure  4-4  illustrates  better  the  nodal  spacing  as  a  function  of  axial  distance  given  at  the 
early  time  range  of  T=0.12  to  T=0.30.  In  this  figure,  you  can  see  a  more  refined  transient 
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behavior.  In  the  preceding  chapter,  it  was  shown  that  the  initial  swelling  of  the  vortex 
tube  occurred  at  approximately  T=0.21,  for  the  current  x=0.20  case.  In  agreement  with 
this,  it  is  clearly  evident  from  Figure  4-4  that  pile-up  has  also  started  at  this  time.  As 
continuity  states,  this  supports  the  argument  that  axial  pile-up  and  radial  expansion  begin 
simultaneous. 


Axial  Node  Pile-Up  Investigation 
Similar  Linear  Velocity  and  Circulation  Ramp-Up 
x  =  0.2  sec  Cmax  =  20  Umax  =  5 


Figure  4-4:  Nodal  spacing  (Az)  vs.  axial  distance,  short  range 
4.1.2  Radial  Pile-up 

The  previous  section  investigated  the  axial  movement  of  the  filaments  subject  to 
induced  velocities  in  relation  to  each  other.  In  doing  such,  the  connection  between  the 
filament  behavior  during  transient  breakdown  and  that  of  a  shockwave  was  made. 
However,  the  results  can  be  deceiving  because  the  investigation  looked  only  at  the  pile- 
up  of  nodes  along  the  streamwise  axis,  z,  excluding  any  pile-up  effect  in  the  radial 
direction.  In  order  to  fully  visualize  and  account  for  three  dimensional  pile-up  effects,  we 
must  further  investigate  the  transient  behavior  of  pile-up  in  the  radial  direction.  One 
method  of  illustrating  the  relative  radial  pile-up  is  to  view  a  meridional  projection  of  one 
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filament,  initially  emerging  at  z=0  at  r=0.03,  the  inner  most  radius,  focusing  on  its 
corresponding  nodes.  An  example  of  this  projection  is  shown  is  Figure  4-5,  which  the 
gives  better  insight  into  the  existence  of  radial  pile-up.  However,  since  it  is  a  meridional 
projection  it  is  still  quite  hard  to  distinguish  between  the  axial  spacing  and  radial  spacing. 
For  example,  from  the  three-dimensional  rendering  representation  in  Figure  4-5,  we 
know  that  the  yellow  filament  traces  a  spiral  streakline  as  it  turns  around  the  vortex  core 
and  into  the  plane,  becoming  hidden  in  this  side-view.  However,  this  fact  is  not  so  clearly 
evident  in  the  meridional  projection  of  the  filament  nodes.  In  fact,  the  two-dimensional 
projection  does  not  aid  in  the  ability  to  distinguish  when  the  filament  turns  into  the 
meridional  plane  from  when  it  turns  out  of  the  plane. 

Vortex  Filament.  Tau :  0-20,  Time  :  1.40 


1 


Figure  4-5:  Radial  pile-up:  meridional  nodal  projection 


This  shortcoming  of  the  meridional  projection  is  evident  throughout  the  entire  history  of 
the  vortex  flow,  as  shown  in  Figure  4-6. 
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Vortex  Filament  Nodes,  Meridional  Plane,  Tau  :  0.20,  Time  :  0.80-2.00 


Figure  4-6:  Radial  pile-up:  time  history  of  meridional  nodal  projection 

Another  method  of  illustrating  the  radial  pile-up  is  to  project  the  filament  nodes 
onto  a  plane  orthogonal  to  the  meridional  plane.  You  can  think  of  this  case  as  an  end- 
view,  sitting  on  the  centerline  far  downstream  and  looking  up  the  vortex  core  (z-axis) 
toward  the  vorticity  source.  This  projected  view  is  illustrated  in  Figure  4-7.  This  view 
gives  a  much  better  interpretation  of  the  radial  pile-up,  as  you  can  now  witness  the 
spiraling  streakline.  In  this  case,  movement  in  the  positive  y-axis  direction  represents  a 
turn  into  the  meridional  projected  plane.  However,  in  opposition  to  the  meridional 
projection,  this  orthogonal  projection  gives  no  information  in  regards  to  the  axial  spacing. 
In  fact,  Figure  4-7  gives  no  information  concerning  the  nodal  overlapping,  which  can  be 
seen  in  the  meridional  plane. 


Vortex  Filaments,  Tau  :  0.20,  Time  :  1.40 


Figure  4-7:  Radial  pile-up:  end-view  projection 
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Although  the  orthogonal  projection  method  offers  additional  spacing  information,  it  is 
no  better  in  representing  the  absolute  spacing  between  nodes  along  the  filament  than  that 
of  the  meridional  projection. 

4.1.3  Pile-up  Along  Filament  Length 

As  shown  in  the  previous  sections,  a  shock-like  pile-up  behavior  occurs  in  both 
the  axial  and  radial  directions.  Therefore,  in  order  to  capture  the  absolute  spacing  of  the 
nodes  we  must  be  able  to  interpret  the  respective  directional  spacing  simultaneously, 
rather  than  independently.  The  axial  projection  and  end-view  projection  alone  fail  to 
illustrate  this  absolute  spacing  between  two  adjacent  nodes,  because  neither  is  able  to 
account  for  simultaneous  movement  in  both  planes. 

Thus,  one  simple  planer  projection  of  node  locations  is  not  sufficient  to  achieve  this 
multi-dimensional  assessment.  Instead,  we  must  represent  the  nodal  spacing  along  the 
vortex  filament  length.  A  perspective  view  is  offered  in  Figure  4-8.  On  the  left  side  of 
the  figure  is  a  three-dimensional  rendering  of  the  r=0.03  streaklines.  Amidst  the  solid 
surface,  a  yellow  line  traces  out  an  individual  vortex  filament.  On  the  right  side,  this 
individual  streakline  is  then  isolated  and  its  respective  filament  nodes  are  highlighted.  (It 
does  not  matter  which  filament  within  the  radial  shell  we  choose  to  highlight,  since  the 
vortex- filament  method  assumes  that  the  spiraling  flow  is  axisymmetric.) 


Vortex  Filaments,  Tau  :  0.20,  Time  :  1.40 


Figure  4-8:  Absolute  pile-up:  perspective  view 
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Figure  4-8  is  able  to  illustrate  the  connection  between  axial  and  radial  pile-up  much 
better  than  any  projected  method.  As  you  can  see,  as  the  axial  pile-up  is  accentuated,  the 
fdament  turns  causing  a  sudden  radial  expansion.  The  filament  then  continues  to  wrap 
around  and  be  drawn  into  the  newly  formed  recirculation  bubble. 

Using  the  nodal  position  data  from  the  vortex-filament  method,  a  plot  of  absolute 
spacing  along  the  vortex  filament  can  be  computed.  This  plot  is  shown  in  Figure  4-9. 
Similar  to  the  axial  pile-up  plots,  the  pile-up  along  the  filament  length  also  exhibits  a 
spatial  compression  and  sharp  singular- like  point,  at  L=3.3.  As  a  means  to  aid  in  relating 
key  positions  along  the  vortex  filament  to  their  respective  axial  locations,  equivalent  axial 
positions  are  noted  in  the  text  boxes.  Since  we  are  investigating  the  spacing  along  the 
filament  length,  it  is  useful  to  note  where  the  filament  obtains  the  farthest  downstream 
position,  prior  to  being  sucked  back  into  the  bubble.  This  position  is  marked  with  a  text 
box  (z=3.45).  Therefore,  any  nodes  further  along  the  filament  length  have  been  sucked 
back  into  the  bubble.  Additionally,  a  side  view  of  the  rendered  vortex  surface  is  also 
shown  in  the  figure  to  help  correlate  regions  of  radial  expansion  with  the  nodal  spacing 
along  the  filament. 


AL  Verses  Vortex  Filament  Length 
t  =  0.20  Time:  1.40  radius  =  0.03 


Figure  4-9:  Nodal  spacing  along  vortex  filament,  r=0.03,  T=1.40 
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As  you  can  see,  the  most  intense  area  of  absolute  pile-up  is  located  at  an  equivalent 
axial  distance  of  3. 14.  This  location  is  the  exact  location  where  the  sudden  radial 
explosion  occurs  at  the  forefront  of  the  recirculation  bubble.  Figure  4-9  only  illustrates 
the  nodal  spacing  along  the  vortex  filament  at  a  particular  moment  in  time  (T=1.40). 
Figure  4-10,  however,  shows  the  temporal  evolution  of  the  nodal  spacing  along  the 
filament.  From  this  figure,  you  can  see  that  nodal  pile-up  begins  early  in  the  transient 
stages  of  breakdown,  where  as  the  sharp  singular- like  point  evolves  at  a  later  time.  This 
sharp  kink  in  the  line  occurs  at  approximately  T=1.20,  which  is  roughly  the  same  time  the 
recirculation  bubble  has  been  fully  established. 


AL  Verses  Vortex  Filament  Length 
t  =  0.20  At  =  0.01  sec  radius  =  0.03 


Figure  4-10:  Time  history  of  nodal  spacing  along  vortex  filament,  r=0.03 
4.2  Vorticity  Gradient  Scenario 

Traditional  views  of  the  vortex  breakdown  have  attributed  the  initial  radial 
enlargement  to  the  presence  of  an  adverse  pressure  gradient,  such  as  the  gradient 
observed  in  a  diverging  pipe  or  over  a  delta  wing.  This  investigation  will  argue  that  the 
initial  vortex  expansion  and  its  subsequent  radial  growth  are  both  a  result  of  the  self- 
induction  among  the  vortex  filaments  subject  to  a  vorticity  gradient.  Furthermore, 
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whenever  the  azimuthal  vorticity  decreases  downstream,  the  resulting  negative 
vorticity  gradient  triggers  both  the  initial  expansion  and  sustained  radial  growth. 


4.2.1  Generation  of  Negative  Azimuthal  Vorticity  Gradient 

Consider  a  vorticity  vector,  <*0 ,  with  end  coordinates  a  and  b  on  the  cylindrical 
surface  of  the  vortex  core  in  the  region  where  the  azimuthal  velocity,  ue ,  is  decreasing 

downstream.  Please  refer  to  Figure  4-11  A  &  C.  Again  we  must  point  out  that  the 
vorticity  vector  always  moves  with  the  fluid.  Therefore,  as  the  vorticity  is  convected 
downstream,  the  downstream  and  trailing  end,  b,  spins  at  a  slower  rate  than  the  upstream 
and  leading  end,  a,  due  to  the  difference  in  the  azimuthal  velocity.  Consequently,  as  the 
vorticity  is  pushed  downstream,  it  transforms  to  <x)' .  This  new  vector,  illustrated  in 
Figure  4-1 1  A,  is  tilted  and  now  defined  with  endpoints  a  ’  and  b  ’.  The  differential  rate  of 
spin  is  the  cause  of  the  titling  behavior,  which  further  causes  the  azimuthal  vorticity 
component  of  <K'  to  be  less  than  that  of  <K . 

Correspondingly,  the  azimuthal  vorticity,  ooe ,  also  decreases  with  downstream 


distance  in  the  region  of  decreasing  azimuthal  velocity,  ue .  The  tilting  behavior 
continues  to  amplify  as  the  vorticity  vector  is  convected  further  downstream.  Eventually, 
the  magnitude  of  the  tilting  gets  so  large  that  ooe  may  become  negative  and  switch 


directions.  Please  refer  to  Figure  4-1  IB.  This  azimuthal  vorticity  sign  switch  may  occur 
even  for  a  vortex  core  of  constant  radius,  somewhere  in  the  region  of  circulation  change, 
or  circulation  ramp-up  in  the  case  of  our  simulation.  In  other  words,  somewhere  between 
the  start  of  the  simulation,  represented  by  point  p,  and  the  time  when  the  maximum 
circulation  value  is  attained,  point  q,  there  exist  a  region  where  the  azimuthal  vorticity 
gradient  becomes  negative. 


5coe 

dz 


<  0  or 


5ooe 

5s 


<0 


The  magnitude  of  negative  azimuthal  vorticity  gradient  in  the  region  of  negative  gradient, 
increases  at  the  same  time  the  sign-switch  occurs.  Please  refer  to  Figure  4-1  ID. 
Accordingly,  this  increase  in  magnitude  further  incites  the  radial  enlargement  of  the 


61 


vortex  core.  It  is  the  inviscid  self-induction  process  subjected  to  a  prescribed  negative 
vorticity  gradient  which  causes  the  initial  expansion  of  the  vortex  tube.  Furthermore,  this 
region  of  initial  expansion,  which  has  been  found  to  occur  in  the  region  of  changing 
circulation,  will  be  shown  to  shift  upstream  into  the  steady-state  region  of  circulation 
with  time. 


Figure  4-11:  Elements  of  vorticity  gradient  generation 
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4. 2. 2  Movement  of  Negative  Azimuthal  Vorticity  Gradient  Region 

As  the  region  of  initial  radial  expansion  further  expands,  it  extends  upstream. 

This  expanding  region  will  continue  to  extend  upstream  and,  eventually,  will  reach  the 
point  of  steady  circulation.  At  this  location,  the  extending  expansion  region  will  further 
induce  a  negative  azimuthal  vorticity  gradient.  This  behavior  is  shown  in  Figure  4-12. 

As  the  fluid  particle  at  a  radius  r  is  pushed  radially  outwards  to  r  +  8r ,  the  downstream 
end  of  (K ,  b,  is  also  radially  displaced.  At  this  point,  the  circulation  is  still  a  constant 
Cmax  value  and  therefore,  the  azimuthal  velocity,  ue ,  at  point  b  is  less  than  that  at  point  a. 

This  further  causes  a  differential  rate  of  spin,  similar  to  the  situation  discussed  in  section 
4.2.1.  Asa  result  of  this  differential  rate  of  spin,  the  azimuthal  component  of  the 
vorticity  vector,  at  the  next  time  step,  is  decreased: 


Figure  4-12:  Elements  of  vorticity  gradient  movement 


Relating  back  to  the  previous  section,  this  decrease  in  the  azimuthal  component  of 
vorticity,  causes  the  negative  vorticity  gradient  to  shift  further  upstream  toward  the 
constant  state  region.  This  shift  triggers  the  newly  occupied  region  to  radially  expand, 
which  then  influences  further  upstream  regions.  Therefore  as  you  can  see,  the  feedback 
cycle  continues  and  the  negative  vorticity  gradient  region  repeatedly  shifts  upstream. 
Here  by  upstream  shift,  it  is  meant  that  the  movement  is  relative  to  the  downstream  end 


63 

of  the  vortex.  Since  the  vortex  is  convected  downstream  in  reference  to  the  absolute 
coordinate,  the  movement  of  the  negative  vorticity  gradient  region  may  be  still  in  the 
downstream  direction. 

4.2. 3  Supporting  Simulation  Results 

The  foregoing  scenario  indicates  that  radial  expansion  occurs  in  the  region  where 
the  azimuthal  vorticity  gradient  is  negative.  Furthermore,  in  a  coordinate  system  moving 
the  with  the  first  vortex  nodes,  this  negative  vorticity  gradient,  and  hence,  region  of  radial 
expansion,  continually  shift  upstream  as  time  progresses.  However,  in  order  to  justify 
this  scenario,  one  needs  more  detailed  numerical  evidence.  We  can  offer  supportive 
computational  evidence  of  the  aforementioned  claims  with  the  use  of  the  vortex-filament 
simulation.  Using  the  results  from  the  simulation,  we  are  able  to  calculate  and  relate  the 
values  of  the  azimuthal  vorticity,  radial  velocity,  and  azimuthal  vorticity  gradient  as  a 
function  of  axial  distance.  This  information,  gathered  for  T=0.21,  is  shown  in  Figures  4- 
13  and  4-14. 


d(DQ/dz  and  co0  Verses  Axial  Distance  At  Constant  Radius 
t  =  0.2  Time  =  0.21  Constant  Radius  =  0.03 


Figure  4-13:  Azimuthal  vorticity  gradient  at  constant  radius,  T=0.21 


d(oQ/dz 


doo0/fe  and  coe  Verses  Axial  Distance  Along  Projected  Streamline 
t  =  0.2  Time  =  0.21  R=0.03 
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Figure  4-14:  Azimuthal  vorticity  gradient  along  streamline,  T=0.21 

Since  we  are  interested  in  the  eventual  control  of  vortex  breakdown  at  the  early  transitive 
stages,  the  early  time  of  T=0.21  was  selected.  Additionally,  this  time  is  also  noted  as  the 
time  when  initial  radial  expansion  is  observed  inside  the  vortex  core.  The  difference 
between  the  two  figures  is  that  values  for  Figure  4-14  were  calculated  for  a  constant 
radial  position  of  r=0.03,  projected  on  the  meridional  plane  ,  where  as  in  Figure  4-15,  the 
values  were  calculated  along  the  projected  streamline  initially  emanating  from  r=0.03. 
Regardless,  both  figures  show  a  positive  radial  velocity  for  the  majority  of  the  negative 
vorticity  gradient  region.  In  fact,  the  only  slight  difference  occurs  at  the  upstream  end  of 
the  negative  vorticity  region.  At  this  location,  there  is  a  slight  delay  in  the  radial 
expansion  after  the  vorticity  gradient  obtains  a  negative  value.  This  lag  is  not  only  small, 
but  expected  when  considering  the  numerically  approximation  of  a  two-phase  flow.  It  is 
also  important  to  note  that  the  maximum  radial  velocity  occurs  when  the  negative 
vorticity  gradient  reaches  a  maximum  negative  value.  The  maximum  occurs  at 
approximately  z=0.35,  which  will  be  shown  in  the  following  to  be  located  approximately 
in  the  middle  of  the  red  section 

Similarly,  we  can  use  the  simulation  results  to  illustrate  the  continual  upstream 
shift  of  both  the  negative  vorticity  gradient  and  the  resulting  radially  expanding  region. 
This  upstream  shift  is  illustrated  in  Figures  4-15  and  4-16.  In  both  figures,  the  red  region 
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is  the  region  of  changing  circulation  and  the  green  region  is  the  region  of  constant 
circulation.  As  you  can  see  in  Figure  4-16,  the  initial  expansion  at  the  outer  edge  of  the 
vortex  core,  r=0.42,  does  not  occur  as  early  in  the  region  of  changing  circulation  and  in 
time,  as  does  the  inner  radius,  shown  in  Figure  4-15.  This  is  an  effect  of  the  expansion 
lag.  The  very  onset  of  radial  expansion  occurs  at  the  inner  most  radius  and  then 
propagates  radially  outwards.  It  is  for  this  reason,  that  the  r=0.42  radial  shell  does  not 
experience  any  radial  expansion  until  a  short  time  after  that  of  the  inner  most  radius. 


Figure  4-15:  Movement  of  initial  expansion,  r=0.03 


Time:  0.30 


Time:  0.36 


Figure  4-16:  Movement  of  initial  expansion,  r=0.42 

In  the  figures,  it  is  also  observed  that  as  the  radial  expansion  region  moves  from  red 
region  to  green,  or  upstream  in  reference  to  the  downstream  end  of  the  red  region,  so  does 
the  increasing  curvature  of  the  vortex  streaklines,  outlining  the  radial  tube.  This 
increasing  bending  represents  the  upstream  movement  of  the  negative  vorticity  gradient. 
As  previously  mentioned,  the  gradient  and  expansion  region  feed  on  each  other. 
Therefore,  as  the  radial  expansion  increases,  moving  upstream,  so  does  the  vorticity 
gradient  and  vice  versa. 
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5  Quantitative  Investigation 

In  this  section,  an  equation  displaying  a  simple  dynamical  relationship  between 
the  azimuthal  vorticity  gradient  and  curvature  of  the  projected  instantaneous  streamline 
will  presented  in  equation  form.  A  method  to  numerically  evaluate  the  individual  terms 
in  the  equation  will  be  then  discussed.  Finally,  using  the  results  from  the  vortex- filament 
simulation,  the  relationship  will  be  validated  computationally. 


5.1  Dynamic  Relationship  in  Meridional  Plane 

In  the  preceding  chapter,  it  was  shown  that  the  initial  expansion  and  subsequent 
growth  of  the  vortex  tube  was  induced  by  the  global  presence  of  a  negative  azimuthal 
vorticity  gradient.  As  a  complement  to  that  qualitative  investigation,  we  offer  here  a 
local  and  dynamical  relationship  between  the  azimuthal  vorticity  gradient  at  a  particular 
location  and  the  curvature  of  the  instantaneous  streamline,  projected  onto  the  meridional 
plane,  at  the  same  location.  In  the  meridional  plane,  the  streamline  coordinates  (s,  n)  are 
defined  such  that  s  is  along  and  n  is  normal  to  the  projected  instantaneous  streamsurface. 
For  the  remainder  of  this  report,  this  projected  surface  will  be  referred  to  as  the 
meridional  streamline.  These  coordinates,  along  with  their  associated  angles  are  shown 
in  Figure  5-1.  Applying  the  inviscid  equations  of  motion  to  this  meridional 
coordinate  system,  one  can  show 
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ooe  is  the  azimuthal  vorticity,  t  the  time,  rm  the  radius  of  curvature  of  the  meridional 

streamline,  qm  =  yjxt 2  +  u2  ,  the  meridional  velocity,  oon  =  ^ co2  +  co2  ,  meridional 
vorticity,  C(r,  z)  =  C(n,  s)  =  27iru0 ,  the  circulation  at  the  position  z  and  at  radius  r, 

P  ii2 

H  =  —  h - ,  the  total  head  pressure  and,  Om ,  the  relative  angle  between  oom  and  qm . 
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The  sign  convention  for  rm  is  also  shown  in  Figure  5-1.  Please  refer  to  Appendix  A  for 
the  step-by-step  derivation  of  Equation  (5-1). 

If  we  assume  the  special  case  of  steady-flow,  the  time-rate  derivative  term  can  be 
removed.  Integrating  the  resulting  equation  in  terms  of  s  yields 


C  ac _ l_dH 

0  47i2r2qm  5n  qm  dn  ' 


(5-3) 


This  is  reduced  to  the  following  if  n  is  replaced  with  the  stream  function,  \\f ,  by 
dn  =  ,  the  Bragg-Hawthorne  equation: 

rqm 

c  ac  aH 
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Figure  5-1:  Streamline  components  in  the  projected  meridional  plane 
In  Equation  (5-4),  Brown  and  Lopez  {ref.  10)  noted  that  if  circulation  is  held  constant 

ac 

along  the  streamline,  —  =  0 ,  such  as  in  the  steady-state  case  after  ramp-up,  then 

as 
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=  — .  This  result  corresponds  to  the  case  where  the  meridional  vorticity 
d\\f  dv|/  w 

vector,  oom  is  aligned  with  the  meridional  velocity  vector,  qm ,  or  simply  Om  =  0 . 
Utilizing  this  alignment,  they  obtained  an  explicit  expression  showing  that  the 
competitive  relationship  between  the  two  right-hand-side  terms  in  Equation  (5-3)  could 
result  in  a  negative  coe  value  when  the  streamsurface  diverges.  Brown  and  Lopez  were 

the  first  to  point  out  the  prominent  role  of  the  azimuthal  sign-switch  in  the  region  of 
radial  expansion.  Brown  and  Lopez’s  argument  gives  compelling  insight  in  the  steady- 
state  scenario.  However,  it  fails  to  offer  any  explanation  into  the  commencement  of  the 
radial  expansion  or  creation  of  the  azimuthal  vorticity  gradient,  both  of  which  are 
prevalent  and  materialize  during  the  formative  stages  of  breakdown. 

Returning  to  the  present  unsteady  transient  case,  we  rewrite  Equation  (5-1)  as 


dt  i  vm 


—  +  ac  +aH 


where 


1  d  1  dC2 


8k  ds  \  r  qm  dn 
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aH  - 


8s\qm  8n 


In  this  compact  form,  ac  and  aH  are  equivalent  to  the  latter  two  terms  in  Equation  (5-1) 
respectfully.  In  the  transient  case,  these  two  terms  can  become  very  small  in  general 

(C  )2x  . 

when  the  dimensionless  parameter,  B  =  — - is  small.  In  the  definition  of  this 

XX 

parameter,  Cxmax  refers  to  the  maximum  total  circulation,  rc  is  the  vortex  core  radius,  x 
is  the  circulation/velocity  ramp-up  time,  and  Ux,  is  the  freestream  or  convected  velocity. 
The  value  of  B  becomes  smaller,  for  instance,  as  you  decrease  the  ramp-up  time.  By 
making  the  assumption  that  conditions  exist  such  that  B  is  small,  Equation  (5-5)  can  be 
further  reduced. 
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Equation  (5-6)  is  a  simple  expression,  which  affords  a  direct  relationship  for  the 


(5-6) 


azimuthal  component  of  a  complex  three-dimensional  vorticity  vector  at  a  particular 
location,  in  terms  of  the  streamline  curvature  at  that  same  location.  Equation  (5-6)  infers 
that  the  change  in  the  azimuthal  vorticity  along  the  meridional  streamline  has  the  same 
sign  as  the  temporal  change  in  the  radius  of  curvature,  rm.  Therefore,  once  radial 


expansion  commences  in  the  region  of  negative  azimuthal  vorticity, 


<9co0 

ds 


< 0 , and 


Iml  >  0  >  it  continues  to  expand  such  that  the  meridional  streamline  becomes  more  curved 
with  time,  rm2  <  rml .  This  fact  supports  our  initial  argument  that  the  negative  vorticity 
gradient  not  only  initiates  the  radial  expansion,  but  also,  feeds  its  subsequent  growth.  On 

the  contrary,  in  the  region  of  a  positive  gradient,  — -  <  0  ,  and  rml  <  0 ,  the  streamline 

ds 

continues  to  flatten  and  level  toward  horizontal  with  time.  In  this  later  case,  the 
subsequent  streamline  is  marked  by  a  larger  radius  of  curvature,  rm2  >  rml .  The 
difference  between  continual  flattening  (b),  and  increasing  curvature  (a),  is  illustrated  in 
Figure  5-2. 


^  m2 


t  t  +  St 


(a)  negative  azimuthal  vorticity  gradient 
and  increasing  curvature 


f 


(b)  positive  azimuthal  vorticity  gradient 
and  continual  flattening 


Figure  5-2:  Temporal  change  in  streamline  curvature 
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5.2  Inclusion  of  the  ac  Term 

In  the  previous  section,  we  presented  the  fact  that  if  the  B  parameter  was  small, 
we  could  neglect  the  contributions  of  the  ac  and  aH  terms  in  the  Equation  (5-5). 

However,  since  the  value  of  B  is  not  particularly  small  for  our  baseline  case,  B=15,  we 
will  elect  to  retain  this  term  for  completeness.  In  addition,  by  calculating  the  values  for, 
and  investigating  the  trends  of,  ac ,  we  will  gain  additional  understanding  into  the 


dynamics  of  the  flow.  On  aH ,  even  though  it  is  possible  to  have  a  gradient  of  H  in 

inviscid  flow  from  the  following  equation, 

DH  _  1  8P 
Dt  p  8t 

we  will,  however,  choose  to  neglect  aH .  The  reason  is  that  we  are  describing  the  vortex 
breakdown  in  the  language  of  incompressible  vortex  dynamics,  where  pressure  terms  do 
not  appear  explicitly. 

dC 

The  strong  variation  of  — appearing  within  theac  term  can  be  seen  in  Figure  5-3. 

2n 


In  this  figure,  a  series  of  contour  plots  of  local  circulation  (x=0.20)  is  overlaid  with 
assorted  instantaneous  streamlines  (emanating  initially  at  r=0.03,  r=0.18,  r=0.36). 


Local  Circulation,  Meridional  plane,  Time:  0.3 


z 


Figure  5-3:  Local  circulation  and  projected  streamlines,  x=0.20 
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Local  Circulation,  Meridional  plane,  Time:  0.6 
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Local  Circulation,  Meridional  plane,  Time:  0.9 
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Local  Circulation,  Meridional  plane,  Time:  1.2 
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Figure  5-3  (cont.) 


As  you  can  see  in  the  figure,  there  is  a  substantial  change  in  circulation  normal  to  the 
streamlines  inside  the  vortex  core  region  (r=0).  On  the  contrary,  as  you  move  away  from 
the  core,  the  circulation  change  decreases  rapidly.  This  fact  confirms  that,  similar  to  a 
Rankine  vortex,  the  flow  field  outside  of  the  core  can  be  considered  irrotational  and  the 
circulation  attains  a  constant  value.  However,  near  the  core,  it  is  clearly  evident  that  we 
can  not  neglect  the  change  in  circulation  normal  to  the  streamline,  and  therefore,  can  not 
neglect  the  ac  term  in  Equation  5-5. 
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If  we  are  to  numerically  calculate  the  value  of  ac ,  we  need  to  first  express  the 


— term,  within  the  ac  equation,  in  terms  of  the  meridional  trigonometry.  By  doing  so, 
2n 


we  avoid  having  to  numerically  take  a  partial  derivative  of  a  partial  derivative,  which  can 
prove  to  be  both  time  consuming  and  computationally  taxing.  Starting  with  the 
divergence-less  property  of  vorticity  and  using  streamline  coordinate  trigonometry,  one 
can  determine  an  expression  for  this  change  in  local  circulation  normal  to  the  streamline, 


dC 

—  =  2ot(os  .  This  can  be  rewritten  as  Equation  (5-2),  from  which  one  can  derive  the 
dn 


following  expression: 
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The  complete  derivative  of  Equation  (5-7)  is  shown  in  Appendix  B.  Substituting  this 
expression  into  Equation  (5-5)  gives  the  final  dynamic  relationship  equation, 
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5.3  Dynamic  Relationship  Computation 

Up  to  this  point,  we  have  only  derived  the  dynamic  equation  relating  the 
azimuthal  vorticity  and  a  change  in  the  shape  of  a  streamline  at  a  particular  location. 
However,  with  the  use  of  the  vortex-filament  method,  we  can  computationally  validate 
this  relationship  expressed  in  Equation  (5-8)  by  calculating  both  sides  of  the  equation.  In 
order  to  illustrate  the  influences  of  the  two  right  hand  side  terms  in  Equation  (5-8),  we 
will  calculate  their  respective  values  independently. 

5.3.1  Time  rate  of  change  in  radius  of  curvature  (xm) 

The  first  step  in  validating  Equation  (5-8)  is  to  calculate  the  time  rate  of  change  in 
the  radius  of  curvature  at  a  particular  point  along  a  given  streamline.  This  is  achieved  by 
first  projecting  the  meridional  streamline  from  the  three-dimensional  stream  surface.  An 
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example  of  the  projected  streamline  is  represented  by  the  white  lines  in  Figure  5-4 
below,  which  also  shows  azimuthal  vorticity  on  the  meridional  plane. 
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Figure  5-4:  Azimuthal  vorticity  and  streamline  projection,  meridional  plane,  T=1.2 


Once  the  desired  streamline,  at  time  t,  is  obtained,  a  position  along  the  streamline, 
such  as  point  P  in  Figure  5-5,  is  selected  and  a  new  streamline  at  some  later  time,  t  +  8t, 
which  passes  through  P,  is  determined.  In  this  case,  since  we  are  taking  a  time  derivative, 
the  spatial  point  of  interest,  P,  must  remain  unchanged.  We  then  calculate  the  difference 
in  radii  of  curvature  at  P  of  the  two  streamlines,  per  unit  time,  8t.. 


Figure  5-5:  Illustration  of  the  radius  of  curvature  method:  theoretical 

As  can  be  seen  Figure  5-5,  the  streamline  passing  through  point  P  at  t  +  8t  is  not 
necessarily  the  same  streamline  passing  point  Q  at  t  +  8t;  therefore,  the  two  streamlines 
originate  from  different  points,  S  and  S”,  at  the  vorticity  source.  Similar  to  a  schematic 
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sketch  of  Figure  5-5,  Figure  5-6  illustrates  this  streamline  relationship  in  terms  of  our 
simulation  results. 


Meridional  Streamlines,  Tau:  0.20,  Time:  0.29-0.30 


This  process  is  then  performed  for  all  points  along  the  original  streamline.  Therefore,  if 
there  are  n  number  of  points  along  the  original  streamline,  n+1  total  streamlines  must  be 
calculated  to  determine  the  time  rate  of  change  in  radius  of  curvature  along  the  entire 
length  of  the  original  streamline.  It  is  also  extremely  important  to  note  that  our  projected 
meridional  streamline  is  determined  computationally  and  thus,  composed  of  a  finite 
number  of  discretized  points  and  must  be  smoothed  as  best  as  possible  to  avoid  infinite 
radii  of  curvatures  values.  If,  in  fact,  the  streamline  is  not  smoothed  adequately,  a 
particular  point  of  interest,  such  as  P  or  Q  in  Figure  5-6,  may  fall  on  a  perfectly  straight 
line  discretized  segment,  leading  to  an  infinite  radius  of  curvature.  For  this  investigation, 
a  cubic-spline  smoothing  method  was  applied  to  smooth  all  streamlines,  prior  to 
calculating  any  radius  of  curvature.  Additionally,  to  avoid  excessively  large  rm  values, 

such  as  rm  — >  go  corresponding  to  near-straight  line  curvature,  — ,  will  be  used  instead  of 

rm 

rm.  An  illustration  of  the  meridional  streamlines  and  the  solid  rendering  of  the  vortex 
tube  are  shown  in  the  Figure  5-7  (top).  These  same  streamlines,  now  denoted  in  white, 


75 


are  further  shown  superimposed  onto  the  azimuthal  vorticity  field  along  the  meridional 
plane  (middle).  Lastly,  the  bottom  of  Figure  5-7  illustrates  the  plots  of  curvature  verses 
axial  distance  for  both  streamline  (a)  and  (b). 
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Figure  5-7:  Radius  of  curvature  for  projected  streamlines,  t=0.005,  T=0.3 
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5.3.2  ac  term  calculation 

Unlike  the  radius  of  curvature,  the  calculation  of  ac  along  a  specific  streamline, 

at  time  t,  is  much  easier.  In  this  case,  since  we  are  taking  a  spatial  derivative  along 
streamwise  distance,  s,  time  remains  fixed.  As  a  result,  only  one  streamline  is  needed. 
Knowing  all  the  physical  parameters,  we  can  calculate  the  value  of 

u9com  cosOm 
9m 

appearing  in  Equation  (5-7)  for  all  points  on  the  desired  streamline;  and  we  further 
calculate  the  spatial  derivative  along  the  streamline  using  a  center-spaced  differencing 
scheme.  Figure  5-8  and  Figure  5-9  show  the  value  of  ac  as  a  function  of  axial  distance 
for  various  projected  streamlines,  at  T=0.3. 

ac  as  a  Function  of  Axial  Distance  (Meridional  Plane) 
x  =  0.20  At  =  0.01  sec  Time:  0.30  sec 
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Figure  5-8:  ac  as  a  function  of  axial  distance,  x=0.20,  T=0.30 


ac  as  a  Function  of  Axial  Distance  (Meridional  Plane) 
t  =  0.005  At  =  0.01  sec  Time:  0.30  sec 
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Figure  5-9:  ac  as  a  function  of  axial  distance,  t=0.005,  T=0.30 
As  shown  in  both  figures,  the  value  of  ac  along  a  streamline  dramatically  increases  as 
you  move  inside  of  the  vortex  core.  Similarly,  as  shown  in  Figure  5-9,  the  value  of 
ac  for  all  streamlines  also  increases  as  you  decrease  the  circulation  ramp-up  time.  We 

would  expect  this  correlation  between  ac  and  decreased  ramp-up  time,  since  a  faster 
ramp-up  would  create  a  stronger  vorticity  gradient,  initially  at  such  an  early  time  of 
T=0.30.  However,  this  increase  in  ac  with  reduced  ramp-up  time  does  not  mean  thatac 
terms  in  Equation  (5-8)  become  more  predominant  in  comparison  to  the  first  term  on  the 


3 

right-hand  side, - 

at 


n 


Vrmy 


0 

.  In  fact,  as  the  ramp-up  time  decreases,  the  — 

at 


vW 


terms 


become  more  dominant,  as  the  preceding  argument  on  the  B-parameter  suggests;  this 
point  will  be  confirmed  in  5.3.4. 
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5.3.3  — —  term  calculation 

ds 

^CO 

To  calculate  the  — -  term,  we  superimpose  the  meridional  streamline  onto  the 

ds 

meridional  cut  of  azimuthal  vorticity,  coe ,  shown  in  the  plots  of  section  3.2.1.  Examples 
of  the  streamline  projections  and  meridional  cut  of  azimuthal  vorticity  are  shown  in 
Figure  5-4  and  Figure  5-7.  We  then  extrapolate  the  values  of  coe  for  every  point  along 

the  streamline,  and  as  we  did  in  determining  a c ,  using  a  center-spaced  differencing 
scheme  to  calculate  the  spatial  derivative  of  azimuthal  vorticity  along  the  streamline. 

5.3.4  Numerical  results 

Having  calculated  all  the  terms,  the  two  sides  of  Equation  (5-8)  can  be  equated 
and  compared  for  various  ramp-up  times,  x.  A  sample  result  is  illustrated  in  the 
following  graphs.  Figure  5-10,  5-1 1,  and  5-12  illustrate  the  validation  of  Equation  (5-8) 
along  three  different  streamlines  (initially  originating  at  r=0.18,  r=0.36,  and  r=0.60 

c3co 

respectfully)  at  T=0.30  for  x=0.20.  In  the  graphs,  the  left  side  of  Equation  (5-8),— — ,  is 

ds 

represented  by  the  black  dotted  line  averaged  over  the  two  successive  time  interval 
T=0.29  and  T=0.30. 


Average  dco0/ds  Verses  -d/dt  [1/rm]  along  Meridional  Plane 
x  =  0.20  At  =  0.01  sec  Time  Range:  0.29-0.30  sec 
Original  Streamline,  r  =  0.18 
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1.6  1.8  2  0 


Axial  Distance,  z 


Figure  5-10:  Sooe/Ss  verses  axial  distance,  r=0.18,  x=0.20,  T=0.30 


Average  dca0/ds  Verses  -d/dt  [1/rm]  along  Meridional  Plane 
x  =  0.20  At  =  0.01  sec  Time  Range:  0.29-0.30  sec 
Original  Streamline,  r  =  0.36 


Axial  Distance,  z 


Figure  5-11:  dcoe/Ss  verses  axial  distance,  r=0.36,  x=0.20,  T=0.30 


Average  dco0/ds  Verses  -d/dt  [1/rm]  along  Meridional  Plane 
t  =  0.20  At  =  0.01  sec  Time  Range:  0.29-0.30  sec 
Original  Streamline,  r  =  0.60 
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Axial  Distance,  z 


Figure  5-12:  <9oV<9s  verses  axial  distance,  r=0.60,  x=0.20,  T=0.30 


0 

The  green  line  refers  to  the  value  of - 

at 


vW 


only,  where  the  red  line  refers  to  the  sum 


of--® 

at 


Vrmy 


andac  .  By  distinguishing  between  these  two  lines,  as  proposed,  we  are  able 


to  isolate  the  individual  effects  of  two  terms  on  the  right-hand  side.  Additionally,  a 
subplot  of  the  solid  rendering  of  the  respective  vortex  tube  is  offered  in  each  graph  to 

c3co 

associate  the  value  of  — -  at  particular  locations  with  regions  of  radial  expansion  and 

as 

c3co 

contraction.  As  mentioned  earlier,  a  positive  — -  value  indicates  a  region  of  radial 

as 

contraction  and  a  negative  value  indicates  a  region  of  radial  expansion.  Similar  to 

as 


all  previous  solid  renderings,  the  red  region  refers  to  the  region  of  changing  circulation. 

From  the  graphs,  it  is  quite  evident  that  the  results  from  the  vortex-filament 
method  produce  a  good  agreement  between  the  two  sides  of  Equation  (5-8).  As  expected, 
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the  influence  of  ac  is  not  trivial  and  also  shown  to  decrease  with  increasing  radial 
distance.  In  fact,  for  the  r=0.60  streamline  case,  represented  in  Figure  5-12,  the 

Scop, 


d_ 

dt 


Vrm/ 


term  alone  provides  a  good  approximation  to 


Ss 


Although  the  general 


profiles  of  the  red  line  and  dotted  black  line  match  well  in  all  three  graphs,  we  must  be 
more  concerned,  from  a  controls  perspective,  with  the  agreement  in  the  sign  of  the 

c3co 

values,  specifically  at  the  downstream  location  where  — -  =  0 .  This  location  marks  the 

ds 

transition  from  radial  expansion  to  contraction,  which  corresponds  to  the  similar  change 
depicted  by  the  vortex  streaklines  in  the  solid  rendering  subplots.  As  also  observed  in  the 
graphs,  the  agreement  of  this  transition  location,  between  both  sides  of  Equation  (5-8), 
improves  with  the  addition  of  ac  and  as  you  extend  out  radially  from  the  vortex  core.  In 


general,  the  agreement  between  the  values  of — - ,  which  is  determined  using  the  Biot- 

5s 

Savart  law,  and  the  dynamic  relationship  is  quite  favorable,  especially  considering  the 
simplicity  of  Equation  (5-8).  We  must  not,  however,  forget  to  note  the  slight  difference 
in  values  towards  the  downstream  end  of  the  vortex  tube.  This  difference  may  be  the 
result  of  not  accounting  for  the  head  pressure  effects,  contained  in  the  aH  term. 

Figures  5-13,  5-14,  and  5-15  show  the  same  comparison  as  the  previous  three 
graphs  but  for  a  faster  ramp-up,  t=0.005.  All  the  trends  observed  in  the  previous  t=0.20 
case  are  also  seen  in  the  t=0.005  case.  In  fact,  the  additive  effect  of  ac  appears  to 


decrease  faster  with  increasing  radial  distance,  at  the  faster  ramp-up  time.  This  trend  is 
exactly  what  the  B  parameter  suggested  earlier.  In  addition,  the  faster  ramp-up  case 
shows  better  agreement  between  the  two  sides  of  Equation  (5-8),  particularly  at  the 
downstream  end  of  the  vortex  tube.  In  keeping  with  the  head  pressure  argument,  this 
improvement  at  the  downstream  end  would  suggest  that  the  head  pressure  effects 
diminish  with  decreasing  ramp-up  time  value.  However,  before  such  an  argument  can  be 
confirmed,  future  cases  that  include  aH  must  be  investigated. 
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Average  dco0/ds  Verses  -d/dt  [1/rm]  along  Meridional  Plane 
t  =  0.005  At  =  0.01  sec  Time  Range:  0.29-0.30  sec 
Original  Streamline,  r  =  0.18 


Figure  5-13:  dcoe/ds  verses  axial  distance,  r=0.18,  t=0.005,  T=0.30 


Average  dcoe/ds  Verses  -d/dt  [1/rm]  along  Meridional  Plane 
x  =  0.005  At  =  0.01  sec  Time  Range:  0.29-0.30  sec 
Original  Streamline,  r  =  0.36 


Figure  5-14:  dcoe/ds  verses  axial  distance,  r=0.36,  t=0.005,  T=0.30 


Average  5coe/5s  Verses  -did t  [1/rm]  along  Meridional  Plane 
x  =  0.005  At  =  0.01  sec  Time  Range:  0.29-0.30  sec 
Original  Streamline,  r  =  0.60 
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Figure  5-15:  dcoe/ds  verses  axial  distance,  r=0.60,  t=0.005,  T=0.30 

In  summary,  the  vortex-filament  simulation  successfully  validates  the  dynamic 
relationship  expression,  as  well  as  all  findings  originating  from  Equation  (5-8).  Most 
important,  a  simple  method  to  calculate  the  azimuthal  component  of  the  complex  three- 
dimensional  vorticity  vector  in  terms  of  the  temporal  change  in  streamline  curvature  has 
been  identified.  A  general  overview  of  this  relationship  and  its  implications  are  shown  in 
Figure  5-16.  In  this  figure  two  scenarios  are  shown.  The  top  of  Figure  5-16  shows  two 

identical  plots  of  — -  verses  axial  distance  for  x=0.20  at  T=0.30.  On  this  curve  two 


5s 


5oofl 


points  have  been  highlighted:  point  P  with  positive  curvature  belonging  to  -  < 

0  and 

5s 

c3co 

point  Q  with  negative  curvature  belonging  to  -  >0.  Below  both  plots,  the  time 

5s 

evolution  of  both  the  streamline  and  solid  rendering  are  shown  for  T=0.30  to  T=0.40.  As 
you  can  see  on  the  left,  the  radius  of  curvature  at  point  P  along  the  streamline  increases 
from  T=0.30  to  T=0.40,  signaling  a  region  of  radial  expansion.  On  the  contrary,  the 
radius  of  curvature  at  Q,  shown  on  the  right  of  the  figure,  becomes  less  negative 
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( |rmi  |  >  |rm2 1 )  with  time,  signaling  radial  contraction.  The  expansion  and  contraction 


behavior  at  P  and  Q,  respectively,  is  confirmed  by  the  superimposed  solid  renderings. 


Figure  5-16:  Overview  of  dynamic  relationship  expression 
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6  Steadiness  Investigation 

In  this  section  we  will  numerically  trace  the  helical  angles  of  velocity  and 
vorticity  as  the  vortex  tube  transients  downstream.  The  results  will  be  presented  in  a 
temporal  evolution  form  and  will  give  further  insight  into  the  features  of  the  swirling 
vortex  tube.  Additionally,  using  a  steady  state  equation  offered  by  Brown  &  Lopez,  we 
will  attempt  to  assess  the  steadiness  of  our  transient  simulation  by  comparing  the 
vorticity  profile  along  a  respective  streamline  with  that  predicted  by  the  steady  state 
equation 


6.1  Helical  Angles  Calculation 

In  modifying  the  present  vortex- filament  method,  we  determined  for  the  helical 

oc 

angles  of  vorticity  at  its  source  based  on  the  prescribed  —  value  of  1.5.  These  angles  are 


a  function  of  both  radial  distance  and  time:  a  =  — 

u7 


CUfltl-e-1*2) 

2rcrU00(t) 


.  However,  if  the 


velocity  at  any  time  is  proportional  to  that  of  the  circulation,  (^(t)  oc  U00(t) ,  as  is  the 
case  with  our  baseline  calculation,  the  time  dependency  falls  out  and  the  helical  value 
becomes  only  a  function  of  radial  distance.  The  helical  angle  approaches  zero  at  both  the 
center  of  the  vortex  core  and  at  the  far  field  radial  boundary.  These  helical  angles, 
shown  in  Table  6-1,  are  then  applied  to  the  discretized  radial  locations  at  the  vorticity 
source.  As  you  can  see  from  the  values  in  the  below  table,  the  helical  angle  appears  to 
reach  a  maximum  at  the  approximate  edge  of  the  vortex  core  (r=0.30)  and  slowly 
decrease  as  you  extend  out  from  this  edge.  Likewise,  the  helical  angle  also  decreases  as 
you  approach  the  centerline  of  the  vortex  core,  however,  at  a  much  faster  rate.  Although 
we  explicitly  specify  the  initial  helical  angle  at  the  vorticity  source,  the  helical  angle  is 
free  to  develop  on  its  own  as  the  vortex  tube  transients  downstream  from  the  source, 
depending  on  the  velocity  induced  by  the  vortex  filaments  at  every  radial  location  within 
the  vortex  tube. 
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Table  6-1:  Prescribed  helical  angle  at  vorticity  source,  z=0 


Radial  Location 

Helical  Angle,  deg 

r  =  0.03 

29.1 

r  =  0.06 

35.9 

r  =  0.08 

42.5 

r  =  0.12 

49.3 

r  =  0.18 

56.8 

r  =  0.24 

58.8 

r  =  0.30 

59.3 

r  =  0.36 

56.9 

r  =  0.42 

53.6 

The  vortex-filament  method  affords  us  not  only  the  ability  to  physically  trace  the 
helical  angle  of  vorticity,  but  the  ability  to  numerically  determine  both  the  helical  angle 
of  vorticity  and  velocity  as  a  function  of  downstream  location  by  tracking  the  relationship 
between  their  respective  axial  and  azimuthal  components.  With  the  use  of  these  axial  and 
azimuthal  components,  we  are  able  to  view  the  temporal  evolution  of  the  tangent  of  the 
helical  angle  of  velocity,  a ,  and  that  of  the  angle  of  vorticity,  (3 .  This  temporal  evolution 
of  angles,  at  constant  radial  cuts  within  the  vortex  core,  is  shown  in  Figure  6-1. 

Alpha  and  Beta  as  a  Function  of  Axial  Distance 
x  =  0.2  sec  Cmax  =  20  Umax  =  5  Time  =  0.30  sec 
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Figure  6-1 :  Temporal  evolution  of  the  tangent  of  helical  angle  of  velocity  and  vorticity 


alpha,  a  beta,  p  alpha,  a  beta,  p  alpha,  a  beta, 


Alpha  and  Beta  as  a  Function  of  Axial  Distance 
t  =  0.2  sec  Cmax  =  20  Umax=5  Time  =  0.60  sec 
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Alpha  and  Beta  as  a  Function  of  Axial  Distance 
t  =  0.2  sec  Cmax  =  20  Umax=5  Time  =  0.90  sec 


Alpha  and  Beta  as  a  Function  of  Axial  Distance 
x  =  0.2  sec  Cmax  =  20  Umax  =  5  Time  =  1 .20  sec 


Figure  6-1  (cont.) 


Alpha  and  Beta  as  a  Function  of  Axial  Distance 
t  =  0.2  sec  Cmax  =  20  Umax  =  5  Time  =  1 .50  sec 
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Alpha  and  Beta  as  a  Function  of  Axial  Distance 
x  =  0.2  sec  Cmax  =  20  Umax  =  5  Time  =  1.80  sec 


Figure  6-1  (cont.) 

Figure  6-1  offers  a  great  deal  of  information  regarding  the  flow  features  of  the  swirling 
flow.  First  of  all,  the  vorticity  sign-switch  location  which  is  marked  by  the  negative 
value  of  P ,  occurs  prior  to  T=0.30  and  continues  to  move  downstream  with  time. 
However,  it  is  important  to  note  that  the  value  of  a  remains  positive  for  quite  some  time, 
despite  a  negative  p  value,  and  does  not  attain  a  negative  value  until  T=1.80,  well  after 
that  of  the  voriticty  sign-switch.  This  points  to  two  possibilities.  The  first  is  that,  with 
azimuthal  vorticity  always  remaining  positive,  the  axial  velocity  reverses  its  sign  only 
after  T=1 .80,  which  is  inconsistent  with  our  previous  results.  The  second,  and  more 
likely  explanation,  suggests  that  the  reversal  in  axial  velocity  and  reversal  in  azimuthal 
velocity  occur  simultaneously.  If  the  two  flows  reverse  at  the  same  time,  the  value  of  a 
remains  positive. 
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Figure  6-1  also  illustrates  that  the  values  and  over-all  profile  of  a  and  (1 ,  for 
each  respective  radial  cut,  appear  to  match  better  with  time.  In  fact,  after  T=0.60,  the  two 
tangents  of  the  helical  angles  show  striking  agreement  at  each  radii.  This  converging 
trend  is  due  to  Helmhotz’s  law  that  the  vorticity  does,  in  fact,  follow  the  material  line  of 
the  fluid;  thus,  even  the  vorticity  vector  and  velocity  vector  are  not  aligned  initially, 
eventually  they  become  aligned.  In  addition,  the  temporal  evolution  also  shows  that  the 
values  of  a  and  p  reach  a  maximum  value  upstream,  just  aft,  of  the  vorticity  source  and 
slowly  decay  thereafter.  These  maxima  occur  shortly  after  the  inception  of  the  vorticity 
sign-switch.  Once  established,  these  maxima  maintain  a  relatively  constant  value  and, 
unlike  the  vorticity  sign-switch  location,  do  not  move  downstream  with  time.  From  the 
point  of  view  of  controlling  breakdown,  this  time  at  which  these  maximum  values  are 
established  represent  a  critical  stage  in  the  formative  stages  of  breakdown  and  will  play 
an  essential  part  in  preliminary  control  simulations,  to  be  discussed  in  the  following 
chapter. 

In  their  steady  vortex  breakdown  investigation,  Brown  &  Lopez  argued  that  in 
order  for  the  breakdown  to  occur,  a  helix  angle  for  the  velocity  which  exceeds  that  of  the 
vorticity  is  a  required  condition  at  the  early  upstream  location.  Although  our  simulation 
is  a  transient  study  of  breakdown,  this  necessary  flow  condition  is  evident  in  Figure  6-1. 
In  fact,  the  value  of  a  appears  to  maintain  a  larger  value  than  p  up  to  approximately 
T=0.60.  With  this  general  consideration  in  mind,  an  important  question  is  raised;  how 
close  to  a  steady  state  flow  is  our  transient  breakdown  simulation?  Obviously,  we  would 
not  expect  the  initial  flow  leading  up  to  breakdown  to  be  steady,  as  the  flow  features  are 
changing  rapidly,  however,  what  can  be  said  about  the  structure  of  the  flow  once 
breakdown  has  been  achieved? 

6.2  Steady-State  Criterion 

In  order  to  assess  the  steadiness  of  our  transient  breakdown  simulation,  we  will 
compare  the  development  of  the  azimuthal  vorticity  component  along  a  streamline 
obtained  from  our  simulation  with  the  steady  predictions  offered  by  Brown  &  Lopez.  For 
the  case  of  steady,  inviscid,  axisymmetric  swirling  flow,  Brown  &  Lopez  derived  an 
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expression  for  the  azimuthal  component  of  vorticity  in  terms  of  local  upstream 
parameters  and  the  radius  of  the  downstream  stream  surface  only.  This  expression,  for 
both  the  case  of  co9o  =  0  and  a>do  *  0 ,  is  given  below. 
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where  a0  =  ue  /  uz  and  P0  =  coe  /  coz  are  the  tangents  of  the  helical  angle  for  the  velocity 

and  vorticity,  respectfully.  Since  in  our  simulations,  the  azimuthal  vorticity  component 
at  the  vorticity  source  is  defined  by  a  finite  and  non-zero  value,  we  will  focus  on  the  latter 
equation,  where  co9o  ^  0  .  According  to  this  Equation  (6-2),  at  some  location  downstream 

of  a  initial  point,  z0 ,  the  ratio  of  upstream  azimuthal  vorticity  on  the  streamline  to  its 


initial  value  at  z0  is  simply  a  function  of  the  initial  helical  angles  of  the  vorticity  and 
velocity  at  z0  and  the  ratio  of  the  radius  at  the  downstream  location  to  the  initial  point. 

oc 

Brown  &  Lopez’s  claim  that  —  must  be  greater  than  one  for  the  onset  of  breakdown  to 

occur  is  apparent  in  Equation  (6-2).  In  order  for  azimuthal  vorticity  to  switch  its  sign,  a 
required  criterion  for  vortex  breakdown,  corresponding  to  a  diverging  stream 

OC 

surface,  r  >  r0 ,  is  that  the  second  term  in  Equation  (6-2)  must  be  made  positive:  — ^  value 

Po 

greater  than  one.  In  other  words,  the  initial  velocity  helical  angle  must  exceed  that  of  the 
vorticity  in  order  to  spawn  breakdown  since,  from  Equation  (6-2)  as  r  /  r0  increases  from 

a  value  of  one,  only  for  the  case  ofa0  >  P0 ,  co9  /coe  becomes  negative  at  large  r/r0 


values. 
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6.3  Steadiness  Assessment 

Using  Equation  (6-2),  an  initial  assessment  into  the  steadiness  of  our  transient 
vortex  breakdown  simulation  is  conducted.  In  essence,  we  are  comparing  the 
downstream  values  of  the  azimuthal  vorticity,  in  reference  to  some  upstream  reference 
point,  of  our  simulation  results  with  that  offered  by  the  steady  Equation  (6-2).  As 
mentioned,  the  simulation  is  not  expected  to  represent  a  steady  flow  prior  to  the 
formation  of  the  recirculation  region.  This  fact  is  clearly  evident  in  Figure  6-2,  which 
illustrates  the  value  co0  /coe  as  a  function  of  downstream  location  as  calculated  from 

Equation  (6-2)  as  well  as  from  our  transient  simulation  results. 


2.0  T 


Steady  Verses  Unsteady  Investigation 
t  =  0.005  At  =  0.001  sec  Time:  0.18  sec 
Original  Streamline  Y-intersetion  =  0.03 
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Figure  6-2:  Steadiness  assessment:  prior  to  breakdown,  T=0.18 
As  you  can  see,  although  the  two  methods  are  initially  in  good  agreement,  as  they 
should  be  from  the  very  definition  of  coe  /coe  ,  thereafter  coe  /to0  quickly  diverges. 

Furthermore,  it  is  crucial  to  note  that  the  coe  /coe  value  determined  by  Brown  &  Lopez 

always  remains  positive;  this  suggests  that  the  flow  never  experiences  a  vorticity-sign 
switch.  On  the  other  hand,  the  plot  of  coe  /ooe  offered  by  the  transient  simulation  shows 


93 


a  negative  sign-switch  at  approximately  z=0.53.  This  location  matches  well  with  the 
sign-switch  location  of  the  vortex  streaklines  seen  in  the  solid  rendering.  Having 
determined  that  the  simulation  results  leading  up  to  the  formation  of  breakdown  do  not 
represent  that  of  steady  flow,  we  will  now  focus  on  the  flow  at  a  later  time;  specifically 
at  the  time  of  approximate  breakdown,  T=1.0,  and  at  the  last  time  step,  T=2.0.  As 
required  by  Equation  (6-2),  an  initial  upstream  reference  point  must  be  selected  in  order 
to  calculate  the  downstream  development  of  the  azimuthal  vorticity.  An  obvious 
reference  location  would  be  at  the  vorticity  source,  z=0.  However,  in  order  to  minimize 
any  biases  or  negating  effect  caused  by  the  ‘necking’  effect  immediately  after  the 
vorticity  source,  a  reference  location  slightly  downstream  of  the  source  will  be  selected, 
For  the  T=1.0  case,  an  axial  reference  location  of  z=0.5  is  chosen.  We  will  attempt  to 
bridge  the  gap  between  the  vorticity  source  and  reference  location  by  projecting  the 
azimuthal  vorticity  value  upstream  as  well  as  downstream  of  the  reference  location.  This 
location  is  marked  by  the  vertical  dotted  line  in  Figure  6-3  (top). 

However,  as  can  be  seen  in  Figure  6-3  (bottom),  this  reference  location  is  not  the 
same  for  the  T=2.0  case.  Again,  our  simulation  is  a  transient  simulation  and  accordingly 
the  recirculation  bubble  is  continually  convected  downstream.  Therefore,  in  order  to 
keep  a  relative  constant  frame  of  reference  we  need  to  follow  the  same  portion  of  the 
vortex  tube  with  time.  In  the  T=1.0  case,  the  distance  from  the  reference  location  to  the 
end  of  the  recirculation  region  is  approximately  4  units  of  length.  Therefore,  to  maintain 
the  same  vortex  tube  length  of  4  units,  a  reference  location  of  z=2.0  is  selected  for  the 
T=2.0  case.  In  other  words,  we  are  only  concerned  with  how  the  downstream  portion  of 
the  vortex  selected  at  T=1.0  has  changed  at  T=2.0. 
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Figure  6-3:  Steadiness  assessment:  case  examples 
Figure  6-3  also  illustrates  three  projected  meridional  streamlines  (r=0.03,  r=0.18,  and 
r=0.42).  It  is  along  these  three  streamlines  that  the  value  of  coe  /coe  ,  in  Figure  6-3,  is 

calculated  using  both  the  results  of  the  vortex-filament  method  and  Equation  (6-2). 

Figure  6-4  shows  the  comparison  between  the  simulation  results  and  the  steady 
state  Equation  (6-2)  at  the  approximate  time  of  the  vortex  experiences  breakdown,  T=1.0 
for  the  baseline  parameters  given  in  Table  6- 1 .  As  can  be  seen  in  the  comparison  plots, 
the  two  lines  initially  start  out  in  good  agreement  but  soon  then  after  begin  to  differ.  Of 
the  three  respective  streamlines,  the  r=0.03  streamline  shows  the  best  agreement  while 
the  r-0.42  streamline  shows  the  poorest  agreement. 
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Steady  Verses  Unsteady  Investigation 
x  =  0.2  sec  Cmax  =  20  Umax  =  5  Time  =  1 .00  sec 
r  =  0.03  Streamline 
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Figure  6-4:  Steadiness  assessment  along  individual  streamlines,  T=1.0 
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In  general,  the  difference  between  the  two  becomes  greater  as  you  move  radially  out 
from  the  vortex  core.  In  fact,  for  the  r=0.42  case,  Equation  (6-2)  predicts  the  axial 
location  of  the  vorticity  sign-switch  to  occur  20  percent  further  upstream  than  that  of  our 
simulation  results.  From  this  information,  we  can  make  the  claim  that  the  fluid  at  the 
very  center  of  the  vortex  core  of  our  transient  simulation  carries  more  quasi-steady  flow 
qualities  than  that  at  the  outer  edge  of  the  core.  Additionally,  it  is  also  evident  from 
Figure  6-4  that  the  agreement  in  the  predicted  location  of  the  vorticity  sign-switch  also 
worsens  as  you  move  out  radially  from  the  vortex  core.  Moreover,  the  gross  trends  for 
the  results  of  both  methods  appear  to  be  relatively  similar  along  all  three  streamlines.  In 
fact,  both  methods  illustrate  a  vorticity  sign-switch  along  each  respective  streamline, 
which  we  know,  from  our  previous  studies,  to  occur  in  fat  at  this  time.  The  location  of  the 
vorticity  sign-switch,  predicted  by  the  simulation  results,  remains  relatively  unchanged  in 
axial  location  with  increasing  radii.  This  trend  ties  into  our  earlier  argument  in  that  the 
axial  flow  stagnation  point  first  appears  at  the  inner  most  core  radius  and  prorogates 
almost  instantly  radially  outward. 

Continuing  on  with  our  steadiness  assessment,  we  perform  at  T=2.0  a  similar 
comparison  between  the  development  of  azimuthal  vorticity  along  the  streamline  with 
that  proposed  by  Brown  &  Lopez’s  steady  formula,  which  is  illustrated  in  Figure  6-5. 
Immediately  clear  is  the  fact  that  the  three  plots  appear  to  show  better  agreement  than 
those  offered  for  T=1.0  in  Figure  6-4.  Like  the  T=1.0  case,  the  r=0.03  case  at  T=1.0, 
shows  remarkable  agreement  between  the  two  methods  up  to  the  sign-switch  location 
which,  in  reference  to  Figure  6-4,  coincides  with  the  axial  velocity  stagnation  point. 
However,  immediately  after  this  point,  the  two  methods  begin  to  differ  with  the  Brown  & 
Lopez  method  showing  a  faster  negative  drop  in  the  value  of  coe  /ooe  .  Similarly,  the 

profiles  of  the  intermediate  and  outer  streamline  at  T=2.0  maintain  a  better  agreement  for 
a  longer  portion  of  the  vortex  tube.  It  is  also  important  to  note  that  at  r=0.42,  the  Brown 
&  Lopez  formula  does  not  show  the  value  of  coe  /ooe  becoming  negative;  hence 

according  to  their  result  the  vorticity  along  this  material  streamline  does  not  experience  a 
sign-switch,  in  contrast  to  the  results  of  the  simulation.  In  other  words,  the  steady 
formula  suggests  that  the  flow  at  this  radius  continues  to  swirls  around  the  core  and 
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recirculation  bubble  in  its  initial  direction,  never  experiencing  a  switch  in  swirling 
direction. 


Steady  Verses  Unsteady  Investigation 
x  =  0.2  sec  Cmax  =  20  Umax  =  5  Time  =  2.00  sec 
r  =  0.03  Streamline 
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r  =  0.18  Streamline 
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Figure  6-5:  Steadiness  assessment  along  individual  streamlines,  T=2.0 
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In  summary,  the  steadiness  assessment  of  our  transient  vortex  simulation  has 
illustrated  a  trend  opposite  to  that  of  the  radius  of  curvature  method  offered  in  the 
previous  chapter.  The  radius  of  curvature  method  showed  better  agreement  with  the 
simulation  vorticity  results  with  increases  radial  distance  from  the  core  centerline  and 
with  decreasing  time.  On  the  contrary,  the  vorticity  relationship  represented  by  the 
Brown  &  Lopez’s  steady  formula  was  found  to  be  in  better  agreement  with  our 
simulation  results  with  decreasing  radial  distance  and  increasing  time.  As  a  result,  there 
is  a  trade-off  between  the  two  methods;  as  one  method  improves  with  a  varying  degree, 
the  other  worsens.  From  a  control  prospective,  although  the  Brown  &  Lopez  method  does 
give  support  to  the  flow  features  of  swirling  flow,  it  does  not  offer  much  accurate  insight 
into  the  flow  properties  during  the  formative  stages  of  vortex  breakdown. 
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7  Preliminary  Control 

Having  gathered  considerable  great  insight  into  the  flow  mechanisms  of  swirling 
flow  and  isolated  key  features  of  inviscid  vortex  breakdown,  this  section  investigates  the 
possibility  to  control  and/or  delay  vortex  breakdown  by  manipulating  the  azimuthal 
vorticity  gradient.  Using  the  results  presented  earlier  in  this  paper  along  with  past 
contributions  to  the  explanation  of  vortex  breakdown,  this  section  offers  two  specific 
control  simulation  scenarios.  The  first  involves  adding  opposite  circulation  to  the  flow 
whereas  the  second  scenario  entails  creating  a  positive  vorticity  gradient.  Both  control 
cases  will  be  modeled  using  the  vortex-filament  method  and  the  results  will  be  presented 
and  discussed. 

7.1  Mirrored  Profile  Control  Simulation 

The  first  attempt  at  actively  controlling  vortex  breakdown  entails  the  addition  of 
opposite  circulation  into  the  vortex  core.  Although  this  technique  has  been  proven  to  be 
quite  difficult  to  perform  experimentally,  the  vortex-filament  simulation  can  perform  this 
rather  easily  and  consistently. 

7.1.1  Simulation  profile 

But  when  is  the  appropriate  time  to  add  opposite  circulation  to  be  added? 
Obviously,  we  want  to  give  the  simulation  enough  time  to  induce  the  formative  stages  of 
breakdown  such  that  we  can  see  the  effects  of  adding  opposite  circulation.  One 
distinguishing  feature  that  was  discovered  in  the  previous  chapter  was  the  unique 
maximum  value  of  P .  This  maximum  value  occurred  after  the  vorticity  sign-switch  but 
prior  to  the  complete  development  of  the  recirculation  region.  Therefore,  we  will  use  this 
feature  as  the  ‘control  switch’  for  our  simulation.  In  other  words,  when  the  value  of  p 
(along  the  r=0.18  plane)  reaches  a  maximum,  the  simulation  will  immediately  begin  to 
add  opposite  circulation,  equal  to  Cfinai  in  magnitude,  to  the  vortex.  We  will  refer  to  this 
control  simulation  as  the  ‘mirrored’  control  simulation  since  we  are  mirroring  the 
circulation  profile. 
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Ideally,  you  would  want  to  add  the  opposite  circulation  in  order  to  affect 
directly  the  region  of  flow  where  the  formative  features  of  breakdown  are  occurring.  One 
method  to  achieve  this  would  be  to  allow  the  region  of  opposite  circulation  to  be 
convected  downstream  at  a  faster  rate  than  the  rest  of  the  vortex  tube,  allowing  it  to 
catch-up  with  the  formative  breakdown  region.  The  present  vortex- filament  method 
code,  however,  does  not  incorporate  the  ability  to  do  such.  In  our  simulation  every 
material  point  along  the  vortex  is  convected  downstream  at  an  equal  velocity.  Similarly, 
in  our  present  set-up  circulation  is  prescribed  prior  to  the  start  of  the  simulation  and 
therefore,  we  can  not  arbitrarily  change  the  circulation  of  an  individual  portion  of  the 
vortex  downstream  of  the  vorticity  source.  As  a  result,  the  region  of  opposite  circulation 
addition  will  initially  travel  at  the  same  speed  as  the  region  of  initial  breakdown.  This 
should,  however,  not  hamper  the  effects  of  the  opposite  circulation  so  long  as  it  is  added 
close  to  the  formative  breakdown  region. 

Figure  7-1  shows  the  velocity  and  circulation  profile  for  the  mirrored  control  case. 
It  is  also  important  to  note  that  besides  the  switch  in  circulation,  all  other  parameters  are 
set  to  be  the  same  as  the  baseline  case  value,  previously  shown  in  Table  3-1. 
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Mirrored  Simulation:  Circulation  and  Velocity  Profiles 


Figure  7-1:  Mirrored  control  simulation  profile 

7.1.2  Simulation  results 

The  simulation  results  of  the  mirrored  control  case  (r=0.08)  are  shown  below  in 
Figure  7-2.  As  you  can  see,  there  is  a  new  blue  region  associated  with  the  vortex  tube.  In 
an  effort  to  investigate  the  effects  of  radial  expansion  and  contraction  at  the  downstream 
edge  of  the  red  region,  a  vortex  tube  extension  with  no  circulation  was  placed  in  the 
domain  prior  to  T=0.  This  vortex  tube,  represented  by  the  blue  region  below,  is 
analogous  to  the  flow  within  a  tube  pre-existing  prior  to  the  introduction  of  swirling  flow 
or  circulation,  similar  to  the  experimental  transient  study  performed  by  Hamada  {ref.  14) 
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Figure  7-2:  Mirrored  control  simulation  results,  r=0.08 


Mirrored  Profile  Simulation,  Time  :  1.2 
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Figure  7-2  (cont.) 
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7.1.3  Discussion 

The  results  in  Figure  7-2  are  of  arresting  interest.  The  first  thing  to  note  is  that  the 
vortex  tube  undergoes  the  same  exact  formative  breakdown  flow  features  that  were 
discussed  in  the  preceding  chapters.  In  fact,  the  structure  of  the  vortex  tube  is  exact  to 
that  offered  by  the  baseline  simulation  up  to  approximately  T=0.90,  at  which  opposite 
circulation  begins  to  be  added.  The  opposing  circulation  is  evident  by  the  reversal  in 
swirling  direction  of  the  vortex  streaklines,  highlighted  below  for  T=1.0. 


Mirrored  ProWe  StonMion.  T«nw  0.9  MirrOttd  Profile  Si  route)  iOfl.  .  t 


Figure  7-3:  Opposite  circulation  and  vortex  streakline  reversal 
A  sharp  pinch  is  formed  at  the  location  of  the  opposite  circulation  addition 
However,  this  pinch  can  be  explained  from  a  vorticity  gradient  perspective.  The 
circulation  immediately  downstream  of  the  pinch  is  positive  and  the  circulation 
immediately  upstream  is  negative.  As  a  result  of  this  opposite  direction  of  swirl, 
circumferential  tilting  of  vortex  filaments  become  reversed  at  the  kink.  This  generates 
negative  azimuthal  vorticity  upstream  of  the  kink,  followed  by  the  usual  positive 
azimuthal  vorticity  just  downstream.  Consequently,  this  sudden  jump  in  circulation 
causes  a  relatively  large  and  positive  local  azimuthal  vorticity  gradient  at  the  pinch. 
From  our  previous  discussions,  we  know  that  a  positive  azimuthal  vorticity  gradient, 

,  causes  radial  contraction.  However,  in  this  case,  the  gradient  is  so  large  that  the 
dz 

contraction  is  almost  instantaneous,  forming  what  is  illustrated  as  a  sharp  kink  at  T=1.0 
in  Figure  7-2. 
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Subsequently,  another  notable  flow  behavior  becomes  apparent.  At 
approximately  T=1.2,  the  downstream  bowl-like  structure,  which  can  be  thought  of  as  the 
upstream  edge  of  the  recirculation  bubble,  begins  to  lose  its  conical  shape  and  flattens 
vertically.  Quickly  thereafter,  the  once  bowl- structure  flips  and  begins  to  form  upstream. 
Furthermore,  the  addition  of  opposite  circulation  continues  in  time,  the  now  upstream¬ 
facing  structure  begins  to  deepen,  as  seen  at  T=1.40.  However,  if  we  take  a  snapshot  of 
the  flow  just  prior  to  the  flip,  we  see  that  the  flow  exhibits  a  very  thin  radially-shaped 
structure,  as  seen  at  T=  1.30  in  Figure  7-4. 


Figure  7-4:  Mirrored  control  simulation:  vorticity  destruction 
This  disk-like  structure  is  a  result  of  the  vortex-filament  method’s  inability  to  account  for 
viscous  effects.  As  mentioned  earlier,  the  present  simulation  is  inviscid,  and  as  such, 
does  not  allow  for  the  cancellation  of  opposing  vorticity  which,  in  real  flow,  occurs 
within  the  recirculation  region.  As  a  result,  the  induction  process  involves  the  two 
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opposing  vorticity  elements,  inducing  a  velocity  in  the  radial  direction.  Therefore, 
given  an  equal  and  opposite  amount  of  vorticity  within  the  recirculation  bubble,  the 
simulation  would  induce  radial  displacement.  In  the  mirrored  circulation  case,  the 
introduction  of  opposite  circulation  causes  a  mass  inflow  into  the  recirculation  bubble, 
which  from  continuity  and  a  lack  of  viscosity  will  lead  to  larger  radial  displacement. 
This  scenario  is  illustrated  below  in  Figure  7-5. 


Baseline  Case 


No  inflow  due 
to  cancellation 


No  inflow  due 
to  cancellation 


+CO0 


Mirrored  Control  Case 


-OJe  due  to  opposite  circulation 


Figure  7-5:  Vorticity  Cancellation 

The  flip  and  upstream  facing  bubble  is  therefore  the  result  of  the  continuing  induction 
process  on  the  should-be-canceled  region.  These  findings  support  the  claim  that  had  the 
flow  been  given  the  ability  to  cancel  out  opposite  vorticity,  the  radial  expanding  region 
would  have  been  suppressed  and  soon  abolished,  eliminating  any  future  breakdown  in  the 
vortex. 


7.2  Positive  Vorticity  Gradient  Control  Simulation 

The  second  control  scenario  involves  applying  a  positive  vorticity  gradient  to  the 
flow  during  the  formative  stages  of  vortex  breakdown.  The  aim  is  that  such  a  positive 
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vorticity  gradient  will  cause  radial  contraction  along  the  vortex  tube,  suppressing  the 
radially  expanding  region  and  ultimately  preventing  breakdown. 

7.2.1  Simulation  profile 

Similar  to  the  mirrored  control  case,  the  question  of  when  to  apply  the  positive 
vorticity  gradient  again  comes  into  play.  The  placement  of  such  a  gradient  is  crucial.  If 
the  vorticity  gradient  is  applied  too  soon,  then  the  self-induction  process  will  not  trigger 
the  formative  features  of  breakdown,  whereas,  if  the  vorticity  gradient  is  applied  too  late 
(too  far  upstream  of  the  predicted  breakdown  location)  its  effect  on  the  region  of  radial 
expansion  and  reversed  axial  velocity  will  not  be  strong  enough  to  suppress  the  inevitable 
breakdown  of  the  flow.  If  the  vortex-filament  method  could  allow  us,  we  would  choose 
to  convect  the  positive  gradient  faster  downstream  than  that  of  the  downstream  portion  of 
the  vortex  tube,  the  hopes  being  that  the  positive  vorticity  gradient  would  add  to  the 
already  existing  negative  gradient  and  cancel.  This  destructive  cancellation  would 
prevent  radial  expansion  which,  if  left  alone,  will  continue  to  intensify  due  to  the 
feedback  cycle  of  the  self-induction  process,  further  triggering  the  inviscid  breakdown 
progression. 

In  our  case,  since  we  can  not  force  the  positive  vorticity  gradient  faster 
downstream  nor  impose  a  negative  gradient  directly  to  a  downstream  location,  we  are 
hoping  that  the  influence  of  the  upstream  positive  gradient  alone  will  weaken  the 
downstream  negative  vorticity  gradient.  If  this  is  in  fact  the  case,  we  may  not  be  able  to 
completely  eliminate  breakdown,  but  could  suppress  its  development.  In  reality, 
suppressing  breakdown  may  be  as  effective  as  completely  eliminating  it;  if  you  can 
suppress  it  long  enough  the  breakdown  location  is  forced  outside  of  the  area  of  interest. 
This  would  be  analogous  to  delaying  breakdown  over  a  delta  wing  long  enough  such  that 
it  breaks  down  well  aft  of  the  trailing  edge  of  the  wing. 

Since  we  expect  the  effect  of  the  positive  gradient  to  be  less  than  that  of  the 
opposite  circulation,  the  positive  gradient  needs  to  be  added  closer  to  the  moment  of 
initial  radial  expansion.  Therefore,  we  will  set  a  temporal  ‘control  switch’  at  T=0.5. 
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This  time  is  approximately  half  way  between  the  time  of  initial  radial  expansion  and 
the  time  of  maximum  (1  value. 

Figure  7-6  shows  the  velocity  and  circulation  profile  for  this  positive  vorticity 
gradient  control  case. 


Positive  Vorticity  Gradient  Simulation:  Circulation  and  Velocity  Profiles 


Figure  7-6:  Positive  vorticity  gradient  control  simulation  profile 


7.2.2  Simulation  results 

The  simulation  results  of  the  positive  vorticity  gradient  control  case  are  shown  in 
the  following  figures.  Figure  7-7  illustrates  the  solid  rending  results  for  r=0.08,  whereas 
Figure  7-8  shows  the  solid  rending  of  the  vortex  tube  at  r=0.03.  As  in  the  first  control 
simulation,  the  addition  of  a  blue  region  is  utilized  to  model  the  existence  of  fluid  within 
the  domain,  prior  to  introducing  any  type  of  swirling  motion. 


Positive  Gradient  Control  Simulation,  Time  :  0.2 

1  r 

0.5  - 


"  — 

-0.5  - 

-1  - 1 - 1 - 1 - 1 - 1 - 1 - 1 - L 

012345678 


9  10 


Positive  Gradient  Control  Simulation,  Time  :  0.4 


Positive  Gradient  Control  Simulation,  Time  :  0.6 


Positive  Gradient  Control  Simulation,  Time  :  1 


Figure  7-7:  Positive  gradient  control  simulation  results,  r=0.08 
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Positive  Gradient  Control  Simulation,  Time  :  1.2 


Positive  Gradient  Control  Simulation,  Time  :  1.4 


Figure  7-7  (cont.) 


Positive  Gradient  Control  Simulation,  Time  :  0.3 
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Positive  Gradient  Control  Simulation,  Time  :  0.6 


Positive  Gradient  Control  Simulation,  Time  :  0.9 


Positive  Gradient  Control  Simulation,  Time  :  1 .2 


Positive  Gradient  Control  Simulation,  Time  :  1 .5 
1  r 


.'I  _ I _ I _ I _ I _ I _ I _ I _ I _ I _ I 

0123456789  10 

z 


Figure  7-8:  Positive  gradient  control  simulation  results,  r=0.03 
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7.2.3  Discussion 

The  results  in  Figure  7-7  and  Figure  7-8  do  not  illustrate  the  dramatic  deviation  in 
flow  development  as  seen  in  the  previous  control  scenario.  At  first  glance  the  flow 
appears  to  follow  the  exact  formative  breakdown  route  as  our  previous  baseline  scenario. 
However,  through  close  inspection,  we  soon  gain  evidence  that  the  flow  is  different  than 
the  baseline  case.  When  directly  compared  with  the  baseline  case,  there  appears  to  be  a 
slight  delay  in  the  formative  features  of  breakdown  after  the  initial  introduction  of  the 
positive  vorticity  gradient,  commencing  at  T=0.5.  After  this  time,  the  development  of  the 
transient  vortex  is  subject  to  the  effects  of  the  upstream  positive  gradient.  The  initial 
influence  of  the  positive  gradient  on  the  vortex  flow  structure  can  be  seen  in  Figure  7-9. 
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Positive  Gradient  Control  Simulation,  Time  :  0.6 
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Figure  7-9:  Positive  gradient  control  simulation  comparison,  T=0.6 
Figure  7-9  offers  a  direct  comparison  of  the  solid  rendering  of  the  vortex  streaklines  of 
both  the  baseline  and  control  simulation  at  T=0.6.  As  you  can  see,  though  slightly 
different  in  shape,  the  general  profiles  of  the  two  vortex  tubes  agree  well  with  each  other. 
This  strong  similarity  highlights  the  fact  that  the  effect  of  the  positive  gradient  has  been 
quite  minimal  up  to  this  point  in  time.  This,  however,  is  to  be  expected  since  the  positive 
gradient  has  only  been  introduced  for  approximately  0. 1  units  of  time. 
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As  the  flow  continues  to  transient  downstream,  the  vortex  tube  is  subject  to 
more  of  the  positive  vorticity  gradient,  traveling  upstream  of  the  radial  expanding  region. 
Therefore,  we  would  expect  to  see  a  stronger  influence  of  the  gradient  on  the  flow  as  time 
progresses.  Our  assertions  are  proved  to  be  correct  in  Figure  7-10.  In  this  figure,  a 
comparison  between  the  baseline  and  control  simulation  is  again  offered,  but  this  instance 
at  a  later  time  of  T=1.4.  The  structure  of  the  two  expanding  regions  are,  again,  very 
similar;  illustrating  the  initial  formation  of  the  recirculation  bubble.  However,  in  this 
case,  there  is  a  significant  difference  in  the  axial  location  of  the  upstream  edge  of  the 
breakdown  bubble.  The  difference  can  be  traced  back  in  time  to  the  initial  introduction 
of  the  positive  gradient  and  increases  with  time. 
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Figure  7-10:  Positive  gradient  control  simulation  comparison,  T=1.4 
This  behavior  leads  to  the  claim  that  the  addition  of  the  positive  vorticity  gradient  causes 
a  delay  in  the  breakdown  of  the  vortex.  One  question  that  needs  to  be  further 
investigated  is  that  whether  the  positive  gradient  continually  suppresses  the  breakdown 
region  as  it  transients  downstream  or  if,  in  fact,  the  growth  of  the  radial  expansion  region 
is  only  initially  stunted  by  this  gradient  but  not  its  subsequent  growth.  Both  of  these 
cases  could  lead  to  the  flow  development  exhibited  in  Figure  7-9.  Regardless  of  them, 
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the  introduction  of  the  positive  vorticity  gradient  does  appear  to  slow  the  rate  of  radial 
growth,  thus  reducing  the  effects  of  the  reversed  axial  velocity.  Both  aid  in  moving  the 
stagnation  location  further  downstream,  relative  to  the  flow. 

Additionally,  one  can  assume  that  the  positive  gradient’s  ability  to  weaken  the 
radial  expansion  would  only  increase  if  we  were  able  to  either  add  it  directly  to  the  region 
of  radial  expansion  or  force  it  downstream  at  a  faster  rate  than  that  of  the  flow. 
Nonetheless,  the  introduction  of  a  positive  vorticity  gradient  appears  to  delay  the 
breakdown  process,  which  in  many  applications  can  be  viewed  just  as  good  as  preventing 
it 


7.3  Experimental  Control  Considerations 

Both  control  scenarios  offer  valuable  knowledge  into  the  control  features  of  the 
swirling  flow.  However,  in  order  to  validate  the  aforementioned  findings  and  claims,  the 
two  scenarios  need  to  be  confirmed  experimentally.  Doing  such  eliminates  many  effects 
which  are  caused  by  the  limiting  assumptions  of  the  vortex-filament  method. 

Furthermore,  substantiation  of  the  scenarios  experimentally  allows  for  effects  of  viscosity 
to  come  into  play  which,  although  not  required  for  breakdown  to  occur,  is  still  an 
important  feature. 

From  an  experimental  control  perspective  it  may,  however,  be  difficult  to  apply 
either  opposite  circulation  or  a  positive  vorticity  gradient  to  the  downstream  portion  of 
the  vortex.  In  fact,  past  experimental  studies,  such  as  performed  by  Wimer,  have  given 
evidence  that  after  a  certain  distance  the  vortex  core  no  longer  entrains  any  fluid  from 
outside  of  the  core.  In  the  specific  case  of  a  delta  wing,  Wimer  found  this  crucial 
distance  to  be  approximately  10  percent  of  the  chord.  At  this  approximate  distance,  the 
vortex  core  is  completely  formed  and  thereafter  remains  constant  along  the  axial  length 
of  the  vortex  until  reaching  the  location  of  breakdown.  This  fact  serves  to  say  that,  if  left 
undisturbed,  the  elements  needed  for  breakdown  are  established  well  upstream  of  the 
breakdown  location.  In  such  case,  death  of  the  vortex  is  not  only  determined,  but  also 
must  set  into  motion,  near  its  birth. 
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With  this  in  mind,  the  difficultly  in  controlling  the  breakdown  of  the  vortex 
becomes  quite  obvious.  In  order  to  offer  any  possibility  of  success,  one  must  be  able  to 
add  either  opposite  circulation  and/or  impose  a  vorticity  gradient  directly  into  the  core 
during  the  early  development  of  the  vortex  in  order  to  control  its  breakdown.  This  task, 
in  itself,  is  extremely  difficult  since  the  core  size,  as  well  as  the  length  of  the  vortex,  in 
most  experiments  is  extremely  small. 
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8  Current  On-going  Studies 

Even  in  light  of  the  computational  and  theoretical  evidence  gained  by  this  report,  a 
complete  understanding  into  the  fundamental  nature  of  vortex  breakdown  is  still  yet  to  be 
achieved.  In  a  continuing  effort  to  pursue  and  reach  an  all-encompassing  agreement  into 
the  breakdown  of  vortex  flow  and  the  interactions  between  its  key  flow  features, 
additional  investigations  are  currently  on  going.  This  section  offers  two  such  on-going 
investigations.  The  aim  and  background  for  each  investigation  is  briefly  discussed,  in  to 
addition  to  highlighting  their  essential  role  in  future  control  simulations. 

8.1  Necking  Effect 

Although  the  modified  vortex-filament  method  has  been  proven  to  successfully 
model  inviscid  development  of  vortex  breakdown,  there  are  still  many  artifacts  of  the 
current  self-induction  process  which  require  further  refinement  and  improvement.  One 
such  area  for  improvement  involves  the  necking  effect  which  occurs  as  the  vortex 
filaments  emerge  from  the  vorticity  source  into  the  flow  domain.  This  immediate  effect 
causes  the  vortex  tube  to  initially  contract,  giving  the  appearance  of  the  quick  narrowing 
behavior  of  the  vortex  core.  Its  effect  can  also  be  seen  in  the  comparison  plots  of 
Equation  (5-7),  by  the  initial  positive  gradient  immediately  downstream  of  z=0.  The 
necking  effect  is  a  trivial  effect  on  both  the  inviscid  development  of  the  vortex  core  and 
its  effective  breakdown  process.  Nonetheless,  this  effect  represents  an  inaccurate  flow 
feature,  and  does  not  manifest  itself  in  experimental  investigations.  Therefore,  further 
modifications  need  to  be  applied  to  the  current  vortex-filament  method  in  order  to 
eliminate  this  erroneous  flow  feature. 

Before  improvements  can  be  made  to  alleviate  the  necking  at  the  vorticity 
source,  the  cause  for  this  must  first  be  determined.  In  the  current  vortex- filament 
method,  the  geometry  of  the  vortex  tube  is  prescribed  prior  to  the  self-induction  being  set 
off.  The  axial  cross-section  of  vortex  is  composed  of  nine  concentric  circles.  Each  axial 
slice  of  these  nine  concentric  circles  is  the  same,  and  axially  displaced  a  defined  distance, 
dz,  from  the  upstream.  As  a  result  of  this  geometry,  the  vortex  filaments  connecting  the 
nodes  between  the  radial  slices  are  not  only  straight  but  aligned  parallel  to  the  z-axis. 
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Therefore  as  the  filaments  emanate  from  the  vorticity  source,  they  enter  the  flow 
domain  in  this  straight  and  horizontal  configuration.  Once  in  the  flow  domain,  they  are 
subjected  to  the  self-induction  process  and  begin  to  tilt  circumferentially.  They  continue 
to  tilt  more  and  more  until  reaching  a  constant  tilt  angle,  at  which  they  will  stay  til  they 
are  incorporated  into  the  vorticity  sign-switch  region.  It  is  important  to  note  that  the  tilt 
angle,  although  it  remains  temporally  constant  in  the  present  simulation,  is  different  for 
each  respective  radial  shell.  It  so  happens  that  the  necking  effect  occurs  exclusely  in  the 
region  in  which  the  filaments  are  transitioning  from  horizontal  to  the  temporally  constant 
angle. 

If  we  analyze  the  behavior  of  the  vorticity  vector  in  this  region,  it  becomes  evident 
that  the  necking  is  related  to  this  change  in  vorticity  vector.  As  the  vortex  filaments 
begin  to  tilt  from  horizontal,  the  azimuthal  component  of  their  vorticity  vector  also  begins 
to  increase  relative  to  the  upstream  position,  creating  a  positive  azimuthal  vorticity 
gradient.  This  positive  vorticity  gradient,  as  explained  in  the  earlier  chapters,  causes 
local  radial  contraction.  This  radial  contraction  manifests  itself  in  the  form  of  the 
perceived  necking  behavior. 

Therefore,  in  order  to  remove  the  necking  effect,  we  must  eliminate  the  tilt 
transition  of  the  vortex  filaments.  The  most  obvious  way  this  can  be  accomplished  is  by 
forcing  the  filaments  within  the  necking  region  to  stay  at  the  required  tilt  angle.  This 
method,  however,  requires  dramatic  changes  in  the  self-induction  algorithm  which  has 
been  proven  to  be  quite  challenging.  The  alternative,  and  more  appropriate  method, 
entails  prescribing  the  vortex  tube  geometry,  such  that  all  vortex  filaments  along  the 
length  of  the  vortex  tube  are  already  tilted  at  the  needed  angle  prior  to  the  start  of  the  self- 
induction.  This  exact  modification  is  currently  being  implemented  into  the  vortex- 
filament  method. 

A  good  approximation  of  the  required  filament  tilt  angle  is  to  specify  the  helical 
angles,  which  were  discussed  earlier  for  our  baseline  simulation,  earlier  in  this  paper. 
Since  our  simulation  assumes  axi-symmetric  flow,  each  node  within  an  axial  slice  of  the 
vortex  tube  is  simply  rotated  by  an  angle  equal  to  that  of  the  corresponding  helical  angle. 
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Due  to  the  upstream  proximity  of  the  affected  region,  this  modification  is  crucial  to 
future  control  attempts  and  simulations. 

Figure  8-1  illustrates  an  image  of  the  current  modified  vortex  tube.  The  blue 
rings  trace  the  axial  slices  of  the  r=0.42  shell  while  the  red  spiraling  lines  represent  the 
vortex  filaments.  As  mentioned,  the  vortex  filaments  no  longer  run  parallel  to  the 
centerline  of  the  vortex  core,  but  rather  swirl  around  the  vortex  at  the  given  helical  angle. 


Figure  8-1 :  Modified  vortex  tube,  r  =  0.42 


The  upper  illustration  in  Figure  8-2  displays  the  same  vortex  tube,  as  above,  but  this  time 
viewed  from  the  side.  The  lower  illustration  in  the  Figure  8-2,  shows  two  of  the 
concentric  radial  shells  simultaneously;  the  black  filaments  refer  to  the  r=0.24  shell  and 
the  red,  again,  to  r=0.42.  As  can  be  seen,  the  two  respective  vortex  filaments  are  not 
subjected  to  the  same  helical  angle.  In  fact,  the  inner  black  filaments  spiral  tighter  than 
those  of  the  outer  filaments. 


Figure  8-2:  Modified  vortex  tube,  side  view 
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Figure  8-2  (cont.) 


8.2  Additional  Analytical  Study 

Similar  to  the  local  dynamic  relationship  that  was  derived  in  Chapter  5  linking  the 
azimuthal  vorticity  gradient  to  the  local  change  in  radius  of  curvature  along  a  projected 
streamline,  an  expression  connecting  directly  the  radial  velocity  and  the  axial  gradient  of 
azimuthal  vorticity  can  also  be  obtained.  To  obtain  this  direct  relationship,  the  equation 


of  continuity  for  axi-symmetric  flow,  +  =  o  ?  js  substituted  directly  into  the 

r  ft  dz 


azimuthal  vorticity  equation,  oo0  = 


dm  du„ 


dz  dr 


,  resulting  in 


d2ur  1  dm  ur  d2ur  dcofl  /r>  , 

— ^  + L — f  +  — ~  —  — (8-1 

dr2  r  dr  r2  dz2  dz 

In  this  expression  the  axial  gradient  of  azimuthal  vorticity  appears  explicitly  as  the  source 
term  for  radial  velocity.  Therefore  if  there  is  no  azimuthal  vorticity  gradient  present, 
radial  velocity  remains  zero.  This  holds  true  regardless  how  large  the  magnitude  of 
azimuthal  vorticity  may  be.  Once  again,  this  makes  sense  if  we  think  in  terms  of  a  single 
row  of  discrete  vortices.  In  this  case,  a  net  radial  induction  occurs  only  is  the  presence  of 
a  change  in  the  relative  strength  of  vortices,  in  the  horizontal  direction.  If  all  the  discrete 
vortices  carry  the  same  value  of  vorticity,  then  the  radial  induce  velocity  would  be 
cancelled  out  by  the  sum  of  the  neighboring  vortices. 


In  the  absence  of  the  term,  the  left  hand  side  of  Equation  (8-1),  would  be  a 
r 

Laplacian.  In  the  present  case,  the  equation  can  be  solved  for  the  radial  velocity  by  using 
the  Hankel  transform  with  Ji  as  a  kernel  (refer  to  Appendix  C) 
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In  conjunction  with  Equation  (8-4),  it  can  be  concluded  that 

Qi  >0 


From  the  above  derived  relationship,  a  key  conclusion  can  be  made. 
r2  +rt2+(z_z')2 

Noting  that  the  term  Q ,  ( - - - — )  is  always  positive,  Equation  (8-2)  states  that 

2  2rr’ 


when  the  gradient  of  azimuthal  vorticity  is  negative  everywhere  inside  the  vortex  tube, 
there  will  be  a  positive  radial  velocity.  On  the  other  hand,  if  the  gradient  of  azimuthal 
vorticity  is  positive  a  negative  radial  velocity  is  generated,  inducing  radial  contraction.  It 
is  important  to  note  that  the  radial  velocity  at  the  observation  point  (r,  z)  should  be  the 
same  as  determined  using  Biot-Savart’s  law.  In  fact,  the  derivation  of  Biot-Savart’s  law 
and  that  of  Equation  (8-2)  are  intertwined,  such  that  one  can  be  obtained  from  the  other. 
The  powerful  benefit  of  Equation  (8-2),  however,  is  that  its  derivation  is  not  as 
cumbersome  and  is  much  more  straightforward  than  that  of  Biot-Savart’s  law. 
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In  order  to  give  evidence  to  the  aforementioned  relationship,  Equation 
(8-2)  must  be  validated  using  the  results  of  the  modified  vortex- filament  method;  similar 
to  what  was  executed  in  the  computational  validation  of  Equation  (5-1).  Doing  such,  not 
only  confirms  Equation  (8-2)  but,  more  importantly,  lends  strong  support  to  the  many 
claims  presented  throughout  this  vortex  breakdown  investigation. 
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9  Conclusions  and  Recommendations 

9.1  Conclusions 

The  self-induction  theory  of  vortex  breakdown  was  investigated  using  a  refined 
and  modified  version  of  Cain’s  vortex-filament  method.  Through  the  use  its  improved 
simulations,  this  modified  method  has  given  further  evidence  that  vortex  breakdown  is 
initiated  by  the  development  of  an  azimuthal  vorticity  gradient  through  an  inviscid  self- 
induction  feedback  mechanism.  Though  inviscid,  this  self-induction  theory  successfully 
demonstrates  how  the  feedback  mechanism  and  its  central  elements  of  radial  expansion 
pile  up,  and  vorticity  sign-switching  culminate  in  a  turning  point,  the  initial  manifestation 
of  vortex  breakdown. 

Unlike  the  many  past  experiments,  the  simulations  of  the  modified  vortex- 
filament  method,  offered  in  this  report,  have  given  rise  to  and  highlighted  the  paramount 
mechanisms  that  initiate  vortex  breakdown  while  expounding  ways  of  eliminating, 
delaying,  or  weakening  breakdown.  With  this  new  paradigm  of  the  inviscid  self- 
induction  process,  we  have  not  only  taken  big  strides  toward  a  general  and  complete 
understanding  of  breakdown,  but  have  established  a  strong  framework  for  future  methods 
of  controlling  and  exploiting  the  breakdown  of  vortex  flow. 

In  this  robust  breakdown  investigation,  the  many  areas  and  rudiments  of  the 
formation  stages  of  vortex  breakdown  have  been  examined  and  scrutinized.  Although 
some  findings  were  more  far  reaching  and  comprehensive  than  others,  their  combination 
presents  a  more  inclusive  overall  picture  of  vortex  breakdown,  its  many  mechanisms,  and 
how  it  can  be  controlled.  A  summary  of  these  findings  are  highlighted  below. 

1 .  The  modified  vortex  method,  with  its  additive  pitch-up  velocity  profile,  models 
the  explosion  of  the  recirculation  bubble  phenomenon  better  than  any  former 
method,  especially  at  the  inner  most  radius. 

2.  In  the  formative  stages  of  breakdown,  reversed  axial  flow  initiates  at  the  inner 
most  radius.  This  back  flow  leads  to  the  establishment  of  a  downstream 
stagnation  point,  which  strengthens  radially  to  form  the  upstream  edge  of  the 
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recirculation  bubble.  Furthermore  it  is  this  initial  explosion  of  the  inner  core 
which  feeds  the  initial  development  of  the  recirculation  bubble. 

3.  It  is  not  the  magnitude  of  the  azimuthal  vorticity  component  that  determines 
vortex  breakdown,  but  rather  the  existence  of  a  differential  change  in  azimuthal 
vorticity.  In  fact,  given  enough  time,  even  the  smallest  differential  change  will 
trigger  the  inviscid  formative  stages  of  breakdown.  However,  the  stronger  the 
differential  change  (stronger  the  azimuthal  vorticity  gradient)  the  more  vigorous 
the  breakdown  process. 

4.  It  is  the  inviscid  self-induction  process  subjected  to  a  prescribed  negative  vorticity 
gradient  which  causes  the  initial  expansion  of  the  vortex  tube.  Furthermore,  this 
region  of  initial  expansion,  which  was  shown  to  start  occur  in  the  region  of 
changing  circulation,  shifts  upstream  into  the  steady-state  region  of  circulation 
with  time.  As  the  region  of  initial  radial  expansion  further  expands  upstream,  so 
does  it  feed  the  upstream  movement  of  the  negative  azimuthal  vorticity  gradient 
region. 

5.  The  vortex-filament  method  has  given  evidence  to  the  claim  that  there  is  a  unique 
relationship  between  the  time  of  initial  vortex  tube  expansion  and  the  time  of 
complete  breakdown  (determined  by  the  existence  of  a  turning  point).  This 
suggests  that  once  the  inviscid  self-induction  mechanism  causes  a  radial 
expansion,  the  subsequent  process  of  breakdown  is  already  predetermined. 

6.  Our  many  vortex  breakdown  simulations  give  support  to  the  claims  that  the  initial 
formative  stages  of  the  vortex  breakdown  phenomenon  behave  similar  to  the 
developing  evolution  of  flow  properties  in  a  materializing  shock-wave. 

7.  This  investigation  offered  a  local  and  dynamical  relationship  between  the 
azimuthal  vorticity  gradient  at  a  particular  location  and  the  curvature  of  the 
instantaneous  streamline,  projected  onto  the  meridional  plane,  at  the  same 
location.  This  remarkable  relationship  affords  simple  method  for  calculating  the 
azimuthal  component  of  the  complex  three-dimensional  vorticity  vector  at  a 
particular  location,  in  terms  of  the  radius  at  that  location. 
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8.  The  dynamical  relationship  was  validated  using  the  results  of  the  vortex- 
filament  method  simulation.  Moreover,  it  showed  that  a  positive  azimuthal 
vorticity  gradient  causes  radial  contraction  of  the  vortex  tube  while  a  negative 
gradient  indicates  a  region  of  radial  expansion. 

9.  The  vortex-filament  method  afforded  the  ability  to  physically  trace  evolution  of 
and  numerically  calculate  the  helical  angles  of  velocity  and  vorticity  as  a  function 
of  downstream  location.  The  time  evolution  showed  that  the  values  of  a  and 

P  reach  a  maximum  value  upstream,  just  aft  of  the  vorticity  source,  and  slowly 
decay  thereafter.  These  maximum  occur  immediately  after  the  conception  of  the 
vorticity  sign-switch.  Once  established,  these  maximums  maintain  a  constant 
value  and  their  respective  locations  do  not  move  downstream  with  time.  This 
insight  offers  a  possible  control  switch  for  controlling  mechanisms. 

10.  A  steadiness  assessment  of  the  transient  vortex- filament  method  was  performed 
by  comparing  the  development  of  the  azimuthal  vorticity  component  along  a 
streamline  with  the  steady  numerical  predictions  offered  by  Brown  &  Lopez. 

This  assessment  illustrated  that  the  results  of  the  vortex- filament  method 
approaches  the  steady-state  solution  with  time.  Additionally,  it  illustrated  that  the 
trends  obtained  by  Brown  &  Lopez’s  expression  for  azimuthal  vorticity  are 
opposite  of  those  obtained  by  our  local  and  dynamical  expression.  The  accuracy 
of  our  expression  improves  with  increasing  radial  distance  and  decreasing  time, 
while  Brown  &  Lopez  improves  with  decreasing  radial  distance  and  increasing 
time. 

11.  Using  results  obtained  earlier  in  this  investigation  along  with  past  contributions, 
two  specific  control  simulations  were  performed;  the  first  adding  opposite 
circulation  and  the  second  adding  a  positive  vorticity  gradient  to  the  flow.  Even 
in  light  of  the  transient  nature  of  our  simulation,  both  control  cases,  when 
compared  to  the  baseline  simulation,  showed  signs  of  either  eliminating  or 
delaying  vortex  breakdown. 
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9.2  Recommendations 

Based  on  the  findings  in  this  study,  several  suggestions  for  future  work  are 
warranted: 

1 .  Effort  into  the  employment  of  the  modified  vortex  tube  should  continue.  Once 
employed,  the  comparison  plots  offered  in  Chapter  5  need  to  be  recalculated  with 
the  hopes  of  eliminating  the  initial  positive  azimuthal  vorticity  gradient  caused  by 
the  necking  effects. 

2.  In  order  to  further  refine  and  improve  the  vortex-filament  method,  the  method 
should  be  updated  such  that  the  vortex  tube  geometry  is  created  at  the  vorticity 
source  as  a  function  of  time.  Doing  such  would  allow  both  the  circulation  and 
velocity  pitch-up  to  both  be  a  function  of  time;  further  eliminating  the  many 
difficulties  apparent  in  the  present  version  of  the  method.  Additionally,  this 
improvement  would  reduce  the  size  of  the  outputted  data  files  while  reducing  the 
computational  time  required  running  the  simulation. 

3.  An  experimental  study,  using  DPIV  perhaps,  needs  to  be  performed  to  further 
validate  the  local  and  dynamical  expression,  presented  earlier  in  this  investigation. 

4.  Further  simulations  into  the  subjecting  either  a  positive  azimuthal  vorticity 
gradient  and/or  opposite  circulation  to  a  vortex  during  the  formative  stages  of 
breakdown  should  be  explored. 

5.  The  vortex-filament  method  should  be  modified  to  allow  the  additive  regions  of 
opposite  circulation  and  positive  vorticity  gradient  to  transient  downstream  at  a 
faster  rate  than  that  of  the  vortex  tube.  This  would  allow  the  controlling  agent  to 
catch  up  with  the  region  of  initial  radial  expansion,  thus  increasing  its  ability  to 
either  suppress  and/or  delay  breakdown. 

6.  Experimental  investigations  into  vortex  breakdown  at  large  Reynolds  numbers 
need  to  be  performed  in  order  to  further  validate  our  simulation  results  and 
conclusions. 

7.  Experimentally  determine  the  contribution  of  the  aH  term  in  Equation  (5-5). 
Perhaps  with  the  use  of  pressure  sensitive  paint  (PSP). 
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Appendix  A  -  Local  and  Dynamical  Relationship  Derivation 


In  this  section  we  derive  the  relationship  between  the  gradient  of  azimuthal  vorticity  and 
the  radius  of  curvature  of  the  instantaneous  streamline  by  writing  the  principal  equations 
in  the  streamline  coordinate  system.  If  we  project  the  three-dimensional  and 
instantaneous  stream  surface  onto  the  meridional  plane,  we  introduce  a  new  streamline 
coordinate  system.  The  streamline  coordinates  (s,n)  are  define  such  that  s  is  the  distance 
along  the  instantaneous  streamline  and  n,  normal  to  it.  Furthermore,  for  any  quantity,  F, 


SF  fiF  fiF  . 

=  — cos0m  +  —  sin0„ 


fir  fin 


ds 


(A- 1) 


fiF  fiF  .  fiF 

—  =  -— sin0m+  —  cos0m 
oz  on  OS 


(A-2) 


Where  0m  is  the  angle  between  a  tangent  to  the  instantaneous  projected  streamline  and 
the  z-axis.  Additionally,  the  following  continuity  equation  for  axisymmetric  flow 


1  fi(rur)  |  fiuz  _Q 
r  dr  dz 

Can  be  rewritten  in  the  streamline  coordinate  system,  as 

fi(rqm)  fi0m 

ds  fin 

where  qm  is  the  meridional  velocity 

qm=VUr+Uz  , 

Next,  we  write  the  equation  of  motion  as 

fill 


fit 


=  gradH  -fix  <5  =  0 

£2 


(A-3) 
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Where  H  is  the  total  head,  H  =  —  +  — ,  and  <x>the  vorticity.  The  azimuthal  component  of 

p  2 


the  equation  of  motion  may  be  expressed  in  the  streamline  coordinates 

fiC  fiC  A 
— +qm — =0 
dt  m  fis 


(A-7) 


Where  C  is  the  location  circulation.  The  radial  and  axial  components  of  (A-6)  are 


129 


dur  dH 
— L  +  —  + 
dt  dr 
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du  SH 
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If  one  multiplies  (A-8)  by  uz  and  (A- 9)  byur ,  and  then  subtracting  one  from  the  other, 
the  following  is  reached 

dQ„  1 


(A-9) 
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Differentiating  (A-9)  by  s  and  combining  the  result  with  (A-3),  you  reach 


dc oe  1  drm  1  d 
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It  is  important  to  note  that  in  deriving  (A- 10),  the  below  relationship  is  used 

=J_ 

Ss  rm 

Where  rm  is  the  radius  of  curvature  of  the  instantaneous  streamline, 
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Finally,  the  radius  of  curvature  is  defined  such  that  positive  curvature  refers  to  the 
streamline  being  concave  downwards,  and  thus,  negative  when  concave  upwards. 
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Appendix  B  -  Deriving  and  Solving  for  ac  Term 

In  this  section  we  derive  the  expression  for  the  ac  term  which  is  used  in  Equation  (4-4). 
Staring  from  the  governing  vorticity  principle  div<K  =  0  and  expanding  can  obtain 

0  =  JJJdivfj)  dv  =  JJ  <*o  •  §  ds  (B-l) 

where  K  is  a  unit  vector  pointing  outwards.  Since  Equation  (B-l)  involves  a  surface  and 
volume  integral  it  is  important  to  first  define  our  volume  and  respective  surface 
components.  The  volume  of  interest  and  surfaces  are  illustrated  in  Figure  B-l. 


Figure  B-l:  Elements  of  stream  surfaces  and  streamlines 
To  evaluate  the  surface  integral  in  Equation  (B-l),  we  break  down  surface  integrals  to 
each  respective  surface,  shown  in  Figure  B-2,  and  obtain  the  following  expression  for 
each  of  the  three  respective  surfaces 


S3 


S  =  Si  +  s2  +  s3 


Figure  B-2:  Surface  integral  components 
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||  <§!  •  1?!  dSj  =  C 

(B-2) 

|!  <S2-£2  ds2=-(c  +  dc) 

(B-3) 

||  do3  •  I13  ds3  =  27irdn  •  oos 

(B-4) 

where  tos  refers  to  the  streamline  component  of  com  ,  the  meridional 
projection  of  (Sand  dC  corresponds  to  a  change  in  stream  surfaces  separated  by  dn  shown 
in 

Figure  B-l. 


Figure  B-3:  Vector  components  and  angles  of  meridional  streamline 
Using  the  relationships  illustrated  in  Figure  B-3,  we  also  know  the  following  component 
expression  for  vorticity  and  velocity,  for  the  meridional  streamline. 
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®s  =®m  C0S®m 

®m=V®z+°>? 

qm=VUz+Ur 


Having  obtained  the  above  vorticity  and  velocity  component  expressions,  we  can 
substitute  into  each  of  the  three  component  surface  integrals.  When  this  is  done  the 
following  is  reached 

C-(C  +  dC)+27trdn-(os  =0 


From  this,  we  obtain  the  expression  for  the  partial  derivative  of  circulation 
with  respect  to  n. 


5C  o 
—  =  27irco. 

an 


We  can  now  substitute  this  derivative  directly  into  the  expression  for  ac,  Equation 
B-10,  which  was  derived  in  Appendix  A,  as 
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Furthermore,  substituting  Equation  (B-5)  into  (B-10),  gives  the  following 
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Next,  by  noting  C  =  27irue , 
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Finally,  we  can  cancel  like  terms  and  obtain  the  simplified  expression  for  ac, 
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d 

a=  — 

c  as 

It  is  important  to  note  that  value  of  3>m  can  be  determined  by  using  the  combination  of 
the  following  expressions 


ue(omcos0m 
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Appendix  C  -  Additional  Analytical  Expression  Derivation 


Figure  C-l:  Vortex  tube  geometry 


For  a  vortex  tube  of  radius  a,  we  calculate  the  radial  velocity  inducted  at  point  P(r,  z) 
induced  by  the  vorticity  at  the  source  point  Q  located  at  radius  r’  and  separated  axially 
from  the  point  P  by  a  distance  z\.  The  vorticity  is  in  the  axial  direction  and  its  magnitude 
ooe  is  constant  throughout  the  vortex  tube.  Outside  of  the  vortex  tube,  which  is  open- 

ended  at  z=0,  the  vorticity  is  either  at  r>a  or  z>0. 

By  applying  Biot-Savart’s  law,  the  radial  velocity  at  P,  which  is 
axisymmetric  is  give  by 
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Therefore 
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On  r 

ur(r,z)  =  — ^flrlr'dr' 
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where 
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Making  the  change  integral  variable,  Iri  may  be  written  further  as 
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I®  and  l[2)  may  be  expressed,  respectively,  as 

I«=E(k) 


(C-8) 


(C-9) 


I(2)  _ . 
irl 


1 


k2[K(k)-E(k)] 

Substituting  Equation  (C-9)  and  (C-10)  into  Equation  (C-6)  yields 
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Finally,  when  we  substitute  the  above  Equation  (C-ll)  into  Equation  (C-4), 
Equation  (7-2)  in  the  main  text  is  reached 
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Appendix  D  -  Modified  Vortex-Filament  Method  Codes 


D.l  Vlinec_rad_spoke_wfj  Main  Code 

// - 

//  This  code  is  a  modified  version  of  CB  Cain's  original  vortex- filament 
//  method  (2001).  It  was  modified  to  incorporate  multiple  features  to 
//  include  variable  velocity  ramp-up  profile  (temporal  ramp-up)  and  azimuthal 
//  vorticity  calculation  at  the  vorticity  source.  The  code  uses  internal 
//  packages  which  are  required  by  Borland  Builder  C++.  However,  the  code  as 
//  been  restructured  to  such  that  it  can  be  complied  and  executed  by  other 
//  programs,  such  as  Microsoft  Visual  C++. 

// 

//  Update  WFJ  Version  9.1 
//  August  2003 

// - 


//  Internal  Call  Packages 

// - 

#include  <vcl.h> 

#include  <fstream.h> 

#include  <istream.h> 

#include  <iostream.h> 

#include  <conio.h> 

#include  <fastmath.h> 

#include  <time.h> 

// - 

#pragma  hdrstop 

//  Initialize  All  Required  Subroutines 

// 

int  readdata  (); 
int  init  (); 

void  cond_out(char*); 

int  out(double**,  double**,  double**,  char*,  int,  int); 

int  out2(double**,  char*,  int,  int); 

void  calc  ind  (); 

void  r  calc  ind  (); 

void  init_ind(); 

void  u_update(); 

void  sym_calc(); 

void  r_sym_calc(); 

void  init_incase_rad_spoke(int, double**,  double**,  double**, double**, double, 
int, int*, int, int, double*, double, double, double, double, double, 
double, double, double**, double**, double**, double**, int, int, int); 

//  Declare  Global  Variables 

// - 

int  in_case,velvor_on,blow_on,n_fil,wall_on,wall_ax_node,wall_circ_node; 
double  t_fmal,U_inf,beta,r_tubeO,r_tubel,r_tube2,r_tube3,r_tube4,r_tube5, 
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r_tube6,r_tube7,r_tube8,r_tube9,r_tubel0,r_tubel  1  ,r_tubel2; 
double  sigma_init,C, alpha, U; 

double  z_on,z_init,x_init,y_init,dz,dt,t_v_out,x_cent,z_start; 

int  tindex,  *on_node,  max_node,out_skip,t_out_start,out_count,x_nodes, 

z_nodes,end_fix,*r_fil,n_rad,r_fil_in,r_n_fil,r_max_node,r_on_fil; 

double  t_real,pi,*core_con,*r_tube; 

int  n_spokes, counter; 

char  *r_node_out,  *r_sigma_out, * gammaout,  *r_gamma_out, *node_out, * sigmaout, 
*tan_omega_out,  *Utheta_out,  *wz_out,  *ommega_theta_out; 
double  u_induced_x,u_induced_y,u_induced_z,  *U_instant; 
double  *Update_C,  *Radius; 

double  **x_n,  **y_n,  **z_n,  *  *  uinducedtotalx,  **u_induced_total_y, 

**u_induced_total_z,  ** sigma; 

double  * *gamma,start_C_inc,stop_C_inc,C_min; 

double  **r_gamma,**r_x_n,  **r_y_n,  **r_z_n,  **r_u_induced_total_x, 

*  *r_u_induced_total_y,  *  *r_u_induced_total_z,  *  *r_sigma, *  *r_core_con, 
**r_seg_length; 

double  **Utheta,  **wz,  **tan_omega,  **ommega_theta; 
int  main(int  argc,  char  **argv){ 


//  Define  Data  Output  Files 

// 

nodeout  =  ”  vlinec.dat"; 
sigmaout  =  ”  sigma_out.dat'’; 
rnodeout  =  "rvlinec.daf ’; 
rsigmaout  =  ”r_sigma_out.datM; 
gamma  out  -  ’’gamma_out.dat”; 
rgammaout  =  '’rgammaout.daf'; 
tanomegaout  = M  tan_omega_out.dat”; 

Uthetaout  =  "Utheta_out.dat'’; 
wz  out  -  "wz_out.dat"; 

ommegathetaout  =  "ommega_theta_out.dat"; 

int  cl  =  clock(); 
read_data();  //read  input  file 

init(); 

//  Initialize  Program  Variables  (above)  and  output  messages  to  screen 

// - 

cout «  "USE  FOR  CASES  THAT  ARE  ALWAYS  SYMMETRIC  ONLY!"  «endl; 

cout «  "T  fmal:  "«t_fmal  «  endl; 

cout «  "Total  iterations:  "«t_fmal/dt «  endl; 

cout «  "Max  nodes:  "  «max_node  «  endl; 

cout «  "U  inf:  "«U_inf  «  endl; 

cout «  "dz:  "«dz  «  endl; 

cout «  "Initial  condition  case:  "«in_case  «  endl; 

//  Initialize  Output  Files.  Note  that  the  program  is  strcuture  such  that 
//  the  output  files  are  continually  written  over 

// - 

ofstream  outfile(node_out,  ios::trunc); 
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if(!outfile)  { 

cout  «"Error  opening  output  file"«endl;} 
outfile.close(); 

ofstream  sigmafile(sigma_out,  ios::trunc); 
if  (Isigmafile)  { 

cout  «"Error  opening  output  file"«endl;} 
sigmafile.close(); 

ofstream  r_outfile(r_node_out,  ios::trunc); 
if  (!r_outfile)  { 

cout  «"Error  opening  output  file"«endl;} 
r_outfile.close(); 

ofstream  r_sigmafile(r_sigma_out,  ios::trunc); 
if  (Irsigmafile)  { 

cout  «"Error  opening  output  file"«endl;} 
r_sigmafile.close(); 

ofstream  gammafile(gamma_out,  ios::trunc); 
if  (Igammafile)  { 

cout  «"Error  opening  output  file"«endl;} 
gammafile.close(); 

ofstream  r_gammafile(r_gamma_out,  ios::trunc); 
if  (Irgammafile)  { 

cout  «"Error  opening  output  file"«endl;} 
r_gammafile.close(); 

ofstream  tan_omegafile(tan_omega_out,  ios::trunc); 
if  (Itan  omegafile)  { 

cout  «”Error  opening  output  file”«endl;} 
tan_omegafile.close(); 

ofstream  Utheta_outfile(Utheta_out,  ios::trunc); 
if  (!Utheta_outfile)  { 

cout  «”Error  opening  output  fileM«endl;} 

Utheta_outfile.close(); 

ofstream  wz_outfile(wz_out,  ios::app); 
if  (!wz_outfile)  { 

cout  «”Error  opening  output  file"«endl;} 
wzoutfile .  close() ; 

ofstream  ommega_theta_outfile(ommega_theta_out,  ios::app); 
if  (lommega  theta  outfile)  { 
cout  «”Error  opening  output  file"«endl;} 
ommegathetaoutfile .  close() ; 

//  Output  Simulation  Time  to  Screen 

// - 

cout «  Mt_index:  "« t  index  «  ”  ”  «  nt_real:  M«t_real «  endl; 


//  Write  Program  Inputs  to  1st  Row  of  Output  File 
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// 

condout(nodeout) ; 
condout(sigmaout) ; 
condout(gammaout) ; 

//  Write  Vortex  Filament  Node  Position  to  File 

// - 

out(x_n,y_n,z_n,node_out,max_node,n_fil);  //main  vortex  tube 
out(r_x_n,r_y_n,r_z_n,r_node_out,r_max_node,r_n_fil) ;  //radial  spokes 

//  Write  Vortex  Filament  Core  Radius  to  File 

// 

out(sigma, sigma, sigma, sigma_out,max_node,n_fil);  //main  vortex  tube 
out(r_sigma,r_sigma,r_sigma,r_sigma_out,r_max_node,r_n_fil) ;  //radial  spokes 

//  Write  Vortex  Filament  Circulation  to  File 

// - 

out(gamma, gamma, gamma, gamma_out,max_node,n_fil);  //main  vortex  tube 
out(r_gamma,r_gamma,r_gamma,r_gamma_out,r_max_node,r_n_fil) ;  //radial  spokes 


//  Set-Up  Velocity  Profile.  Note  that  the  Circulation  Profile  must  also  be 
//  ajusted  in  the  "incaseradspoke"  file  which  sets-up  the  vortex  tube 
//  filament/node  position. 

//  I  realize  that  this  method  of  implementation  is  rather  old-style,  but  it 
//  allow  for  quick  changes  in  profile  and  allows  for  any  possible  profile. 


//  for  (int  m  =  0;  m  <=  1 ;  m++)  { 

//  U_instant[m]  =  0.0; 

//  Update_C[m]  -  0.0; 

//  } 

for  (int  mm  =  0;  mm  <=19;  mm++){ 
U_instant[mm]  =  0.25  *  (mm+1); 
Update_C[mm]  =  (mm+1); 

} 

for  (int  c  =  20;  c  <=  200;  c++){ 
U_instant[c]  =  5.0; 

Update_C[c]  =  20.0; 

} 


for  (int  c  =  0;  c<=19;  c++){ 
Radius[c]  =  (c*0.1)+0.1; 

} 

for  (int  m  =  0;  m  <=  200;  m++)  { 
for  (int  i  =  0;  i  <=  19;  i++){ 


Utheta[m][i]  =  ((Update_C[m]*(1.0-(exp(-1.0*10.9*Radius[i]* 
Radius[i]))))/  (2.0*3. 14 15*Radius[i])); 
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} 

} 

for  (int  m  =  0;  m  <=  200;  m++)  { 

for  (int  i  =  0;  i  <=  19;  i++){ 

tan_omega[m][i]  =  Utheta[m][i]  /  U_instant[m]; 

} 

} 

int  row=201; 
int  col=20; 

for  (int  m  =  0;  m  <=  200;  m++){ 
for  (int  i  =  0;  i  <=  19;  i++){ 

wz[m][i]  =  (Update_C[m]*10.9*(exp(-10.9*pow(Radius[i],2))))/ 

(Pi); 

} 

} 

for  (int  m  =  0;  m  <=  200;  m++){ 
for  (int  i  =  0;  i  <=  19;  i++){ 

ommega_theta[m][i]  -  wz[m][i]*tan_omega[m][i]; 

} 

} 

out2(tan_omega,tan_omega_out,row,col); 
out2(Utheta,Utheta_out,row,col); 
out2(wz,  wz_out,row,col); 

out2(ommega_theta,  ommega_theta_out,row,col) ; 

// - Begin  main  program  loop - 

while  (t  real  <  tfmal-dt)  { 


calc_ind();  //calculate  induced  velocity  on  vortex  filament  nodes 

tindex  =  tindex  +  1 ; 
t  real  -  t_real  +  dt; 
counter  =  counter  +1 ; 


//write  current  simulation  time  to  screen 
int  c2  =  clock(); 

cout «  "t_index:  "«  t_index«"  t  real:  "«t_real  «  "  run  time:  "  « 
double(c2-c  1  )/CLK_TCK  «  endl; 

//write  output  to  file 

if  (t  index  >=  t_out_start  &&  t  index  ==  out_count*out_skip){ 
cout  «endl; 

cout«MWriting  node  position  to  file,  t_index=  "« t_index«endl; 
cout  «endl; 

cout  «"Code  Updated!  "«endl; 
cout  «endl; 
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cout  «”U_instant  =”«U_instant[counter]«endl; 
cout  «"Final  WFJ  Code...  Updated!  ”«endl; 
cout  «endl; 

//write  vortex  filament  node  position  to  file 
out(x_n,y_n,z_n,node_out,max_node,n_fil);  //main  vortex  tube 
out(r_x_n,r_y_n,r_z_n,r_node_out,r_max_node,r_n_fil) ;  //radial  spokes 

//write  vortex  filament  core  radius  to  file 

out(sigma, sigma, sigma, sigma_out,max_node,n_fil);  //main  vortex  tube 
out(r_sigma,r_sigma,r_sigma,r_sigma_out,r_max_node,r_n_fil) ;  //radial  spokes 

//write  vortex  filament  circulation  to  file 

out(gamma, gamma, gamma, gamma_out,max_node,n_fil);  //main  vortex  tube 
out(r_gamma,r_gamma,r_gamma,r_gamma_out,r_max_node,r_n_fil) ;  //radial  spokes 

outcount  =  outcount  +  1 ;  //index  for  writing  output  to  file 

} 

} 

int  c3  =  clock(); 
cout «  endl; 

cout  «”Total  run  time:  "«  double(c3-cl)/CLK_TCK  «  endl; 

cout «  endl «  ’’Press  any  key  to  continue...”; 

getch(); 

outfile.close(); 

return  0; 

} 

// - 


// - Function  readdata - 

//Read  data  from  vlinec.inp.  Input  file  requires  ’$’  characters  and  end  of  line 

// 

int  read  data  () 

{ 

ifstream  infile(”vlinec_rad_spoke.inp’’,  ios::in); 

if  (!  infile)  return  0; 
infile.ignore(  10000,  '$’); 
infile  » in_case;infile.ignore(1000,  ’\n’); 
infile  » t_fmal;infile.ignore(1000,  ’\n’); 
infile  »  U_inf;infile.ignore(1000,  ’\n’); 
infile  »  dz;infile.ignore(1000,  ’\n’); 
infile  »  dt;infile.ignore(1000,  ’\n’); 
infile  »  max_node;infile.ignore(1000,  ’\n’); 
infile  »  r_fil_in;infile.ignore(1000,  ’\n’); 
infile  »  n_rad;infile.ignore(1000,  ’\n’); 
infile  » t_out_start;infile.ignore(1000,  ’\n’); 
infile  »  out_skip;infile.ignore(1000,  ’\n’); 
infile  »  r_tube0;infile.ignore(1000,  ’\n’); 
infile  »  r_tubel;infile.ignore(1000,  ’\n’); 
infile  »  r_tube2;infile.ignore(1000,  ’\n’); 
infile  »  r_tube3;infile.ignore(1000,  ’\n’); 
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infile  »r  tube4;infile.ignore(1000,  '\n'); 
infile  »  r_tube5;infile.ignore(1000,  ’\n’); 
infile  »  r_tube6;infile.ignore(1000,  ’\n'); 
infile  »  r_tube7;infile.ignore(1000,  ’\n'); 
infile  »  r_tube8;infile.ignore(1000,  ’\n'); 
infile  »  r_tube9;infile.ignore(1000,  ’\n'); 
infile  »  r_tubel0;infile.ignore(1000,  ’\n'); 
infile  »  rtubel  l;infile.ignore(1000,  ’\n'); 
infile  »  r_tubel2;infile.ignore(1000,  ’\n’); 
infile  »  sigma  init;infile.ignore(  1 000,  ,\n'); 
infile  »  C;infile.ignore(1000,  '\n'); 
infile  »  start  C  inc ; infile .  ignore(  1000,  ’\n'); 
infile  »  stop_C_inc;infile.ignore(1000,  '\n'); 
infile  »  C_min;infile.ignore(1000,  ’\n'); 
infile  »  alpha;infile.ignore(1000,  '\n'); 
infile  »  z_on;infile.ignore(1000,  ’\n'); 
infile  »  z_init;infile.ignore(1000,  '\n'); 
infile  »  x_init;infile.ignore(1000,  '\n'); 
infile  »  y_init;infile.ignore(1000,  '\n'); 
infile. close(); 
return  0; 


// - Function  init - 

//Initialize  variables 

// - 

int  init  () 

{ 

r  fil  =  new  int  [n_rad]; 
r  tube  =  new  double  [n  rad]; 

U  instant  =  new  double  [t_fmal/dt]; 

Update_C  =  new  double  [t  final/dt] ; 

Radius  =  new  double  [20]; 
nfil  =  0; 

//initialize  number  of  filaments  in  a  given  radius 

r_fil[0]  =  rfilin; 

r_fll[l]  =  rfilin; 

r_fil[2]  =  rfilin; 

r_fil[3]  =  rfilin; 

r_fil[4]  =  rfilin; 

r_fil[5]  =  rfilin; 

r_fil[6]  =  rfilin; 

r_fil[7]  =  rfilin; 

r_fil[8]  =  rfilin; 

//r_fil[9]  =  r_fil_in; 

//r_fil[10]=  r_fil_in; 

//r_fil[l  1]=  rfilin; 

//r_fil[12]=  r_fil_in; 

//Initiallize  vortex  filament  radii 
rtube  [0]=r_tube0 ; 
rtube  [  1  ]  =r_tube  1 ; 
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rtube  [2]=r_tube2 ; 
rtube  [3  ]=r_tube3 ; 
r_tube  [4]=r_tube4 ; 
r_tube[5]=r_tube5 ; 
rtube  [6]=r_tube6 ; 
r_tube  [7]=r_tube7 ; 
r Jube  [8]=r_tube8 ; 

//r_tube[9]=r_tube9 ; 

//r_tube[  1 0]=r  tube  10; 

//r_tube[  1 1  ]=r  tube  1 1 ; 

//r_tube[  1 2]=r Jube  1 2 ; 

for  (int  j  =  0;  j  <=  n  rad- 1 ;  j++)  { 
nfil  =  r_fll  [j  ]  +n_fil ; 

} 

nfil  =  n  fil  +  1 ;  //add  centerline  filament 

double  **segjength  =  new  double  *[max_node]; 

for  (int  i  =  0;i<max_node;i++){seg  length[i]=  new  double  [n  fil] ; } 

//initialize  simulation  time 
t_real  -  0.0; 
t  index  =  0; 
outcount  -  1 ; 

pi  =  M_PI;//defme  pi 

zinit  =  zinit  +  dz; 

//Initiallize  dynamic  arrays 
core  con  =  new  double  [n  fil]; 

on  node  =  new  int  [n  fil] ; 

Utheta  =  new  double  *[201]; 

for  (int  i  =  0;i<201;i++){Utheta[i]=  new  double  [20];} 
wz  =  new  double  *[201]; 

for  (int  i  =  0;i<201;i++){wz[i]=  new  double  [20];} 
tan  omega  =  new  double  *[201]; 

for  (int  i  =  0;i<201;i++){tan_omega[i]=  new  double  [20];} 
ommega  theta  =  new  double  *[201]; 

for  (int  i  =  0;i<201;i++){ommega_theta[i]=  new  double  [20];} 
gamma  =  new  double  *[max_node]; 

for  (int  i  =  0;i<max_node;i++) { gamma [i]=  new  double  [n  fil];} 
x  n  =  new  double  *[max_node]; 

for  (int  i  =  0;i<max_node;i++)(x_n[i]=  new  double  [n  fil];} 
y  n  -  new  double  *[max_node]; 

for  (int  i  =  0;i<max_node;i++)(y_n[i]=  new  double  [n  fil];} 
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z_n  =  new  double  *[max_node]; 

for  (int  i  =  O;i<max_node;i++){z_n[i]=  new  double  [n  fil];} 

u  induced  total  x  =  new  double  *[max_node]; 
for  (int  i  =  O;i<max_node;i++){u_induced_total_x[i]= 
new  double  [n  fil];} 

u  induced  total  y  =  new  double  *[max_node]; 
for  (int  i  =  O;i<max_node;i++){u_induced_total_y[i]= 
new  double  [n  fil] ; } 

u  induced  total  z  =  new  double  *[max_node]; 
for  (int  i  =  O;i<max_node;i++){u_induced_total_z[i]= 
new  double  [n  fil];} 

sigma  =  new  double  *[max_node]; 

for  (int  i  =  O;i<max_node;i++){sigma[i]=  new  double  [n  fil];} 

//radial  spokes  arrays  and  locations 

int  m_z  =10; 
nspokes  =16; 
r_n_fil  =  10*n_spokes; 
rmaxnode  =  5/0.1; 

r  gamma  =  new  double  *[r_max_node]; 

for  (int  i  =  O;i<r_max_node;i++){r_gamrna[i]=  new  double  [r  n  fil] ; } 
rxn  =  new  double  *[r_max_node]; 

for  (int  i  =  O;i<r_max_node;i++){r_x_n[i]=  new  double  [r  n  fil];} 
r_y_n  =  new  double  *[r_max_node]; 

for  (int  i  =  O;i<r_max_node;i++){r_y_n[i]=  new  double  [r  n  fil];} 
r_z_n  =  new  double  *[r_max_node]; 

for  (int  i  =  O;i<r_max_node;i++){r_z_n[i]=  new  double  [r  n  fil];} 

r  u  induced  total  x  =  new  double  *[r_max_node]; 
for  (int  i  =  O;i<r_max_node;i++){r_u_induced_total_x[i]= 
new  double  [r_n_fil];} 

r  u  induced  total  y  =  new  double  *[r_max_node]; 
for  (int  i  =  O;i<r_max_node;i++){r_u_induced_total_y[i]= 
new  double  [r  n  fil] ; } 

r_u_induced_total_z  =  new  double  *[max_node]; 

for  (int  i  =  0;i<r_max_node;i++)  (r_u_induced_total_z[i]= 

new  double  [r  n  fil] ; } 

r  sigma  =  new  double  *[r_max_node]; 

for  (int  i  =  O;i<r_max_node;i++){r_sigma[i]=  new  double  [r  n  fil];} 
r  core  con  =  new  double  *[r_max_node]; 

for  (int  i  =  O;i<r_max_node;i++){r_core_con[i]=  new  double  [r  n  fil];} 
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r  seg_length  =  new  double  *[r_max_node]; 

for  (int  i  =  O;i<r_max_node;i++){r_seg_length[i]=  new  double  [r  n  fil] ; } 
//Initialize  x_n,  y_n,  z_n 

//Note:  z_n  is  to  the  right  (axial  direction),  x_n  is  up, 

//y_n  is  out  of  page.  When  plotting,  plot  z_n  to  the  right,  x_n  up, 

//and  y_n  out  of  page. 

//subroutine  in  init  incase  rad  spoke.cpp  that  initializes  vortex  filament 
//node  location 

init_incase_rad_spoke(in_case,x_n,y_n,z_n, gamma, C,max_node,r_fil,n_fil,n_rad, 

r_tube,x_init,y_init,z_init,dz,start_C_inc,stop_C_inc,C_min, 

r_gamma,r_x_n,r_y_n,r_z_n,r_max_node,n_spokes,m_z); 


//Intialize  sigma,  find  core  constant  to  maintain  constant  volume  in 
//vortex  filament  segment.  This  allows  simulation  to  handle  vortex 
//stretching 

for  (int  n  =  0;  n  <=  nfil- 1 ;  n++)  { 

for  (int  m  =  0;  m  <=  max_node-2;  m++){ 

sigma[m][n]  -  sigmainit; 

seg_length[m][n]  =  sqrt((x_n[m+l][n]-x_n[m][n])* 

(x_n[m+ 1  ]  [n]-x_n[m]  [n])+(y_n[m+ 1  ]  [n]-y_n[m]  [n])*(y_n[m+ 1  ]  [n]- 
y_n[m]  [n])+(z_n[m+ 1  ]  [n]-z_n[m]  [n])*(z_n[m+l  ]  [n]-z_n[m]  [n])); 

} 

core_con[n]  =  sigma[l][n]*sigma[l][n]*seg_length[l][n]; 

}  ^  . ------- 

//intialize  sigma  for  radial  spokes 
for  (int  n  =  0;  n  <=  r  n  fil- 1 ;  n++)  { 

for  (int  m  =  0;  m  <=  r_max_node-2;  m++){ 
r_sigma[m][n]  -  sigma  init; 

r_seg_length[m][n]  -  sqrt((r_x_n[m+l][n]-r_x_n[m][n])* 

(rxn  [m+ 1  ]  [n]  -rxn  [m]  [n]  )+(r_y_n  [m+ 1  ]  [n] -ryn  [m]  [n] )  * 
(r_y_n[m+l  ]  [n]-r_y_n[m]  [n])+(r_z_n[m+l  ]  [n]-r_z_n[m]  [n])* 
(rzn  [m+ 1  ]  [n]  -rzn  [m]  [n] )) ; 

r_core_con[m][n]  -  r_sigma[m][n]*r_sigma[m][n]* 
r_seg_length[m]  [n] ; 

} 

} 

return  0; 


} 

// - Function  cond  out - 

//  Write  inputs  to  1  st  row  of  output  file 

// 


void  cond_out(char  *out) 

{ 

ofstream  outfile(out,  ios::app); 

//need  to  fill  34  columns 

//output  simulation  variable  to  top  row  of  output  file 
outfile  «  in  case  «  "  "«  n  fil  «  "  "«max_node«  "  M«dt «  " 
«  dz  «  "  M«  Uinf  «  "  "«  t_out_start «  "  "«  outskip  «  "  " 
«  r_tubel  «  "  "«  r_tube2  «  "  "«  r_tube3  «  "  M  «r_tube4«  " 
«  r_tube5  «  "  "«  sigma  init «  "  "«  C  «  "  "«  alpha  «"  " 
«  z_on  «  "  "«  n  rad  «  "  f,«  r_fil[0]  «  "  M«r_fil[l]«M  ” 

«  r_fil[2]«  ”  ”  «r_fil[3]«  "  M«r_fil[4]«  M  M«r_fil[5]«M  M 
«  r_fil[6]«  "  ”  «r_fil[7]«  ”  M«r_fil[8]«  M  M«r_fil[9]«M  M; 
for  (int  p=29  ;  p<=n_fil+l  ;  p++) {outfile  «  0  «  ” 
outfile«  endl; 

} 

// - Function  out - 

//  Write  x_n,  y_n,  z_n  to  file 

// - 

int  out(double**a,  double* *b,  double* *c,  char*out,  int  m  row,  int  n  col) 

{ 

ofstream  outfile(out,  ios::app); 
if  (!  outfile)  return  0; 

//output  x_n 

for  (int  m  =  0;  m  <=  m  row- 1 ;  m++)  { 
for  (int  n  =  0;  n  <=  n  col- 1 ;  n++)  { 
outfile  .precision(6) ; 
outfile  «  a[m][n]  «  " 

} 

outfile  « t_index  «endl; 

} 

//output  y_n 

for  (int  m  =  0;  m  <=  m  row- 1 ;  m++)  { 
for  (int  n  =  0;  n  <=  n  col- 1 ;  n++)  { 
outfile  .precision(6) ; 
outfile  «  b[m][n]  «  " 

} 

outfile  « t_index  «endl; 

} 

//output  z_n 

for  (int  m  =  0;  m  <=  m  row- 1 ;  m++)  { 
for  (int  n  =  0;  n  <=  n  col- 1 ;  n++)  { 
outfile  .precision(6) ; 
outfile  «c[m][n]«"  M; 

} 

outfile  « t_index  «endl; 

} 

//outfile. close(); 
return  0; 


// - Function  out- 

//  Write  tan  omega  data  file 
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// - 

int  out2(double**a,  char* out,  int  mrow,  int  ncol) 

{ 

ofstream  outfile(out,  ios::app); 
if  (!  outfile)  return  0; 

//output  x_n 

for  (int  m  =  0;  m<=m_row-l;  m++){ 
for  (int  n  =  0;  n  <+  n  col;  n++){ 
outfile.precision(6); 
outfile  «  a[m]  [n]  «  "  M; 

} 

outfile  «endl; 

//  for  (int  m  =  0;  m  <=  m  row- 1 ;  m++)  { 

//  for  (int  n  =  0;  n  <=  n_col-l ;  n++)  { 

//  outfile.precision(6); 

//  outfile  «a[m][n]«"  H; 

//  } 

} 

//outfile.  close(); 
return  0; 


// - Function  calc  ind - 

//  Perform  induced  velocity  calculations 

// 

void  calc_ind() 

{ 

init_ind(); 

//calculate  induced  velocity  on  each  node  from  every  other  vortex 
//filament  segment 
int  n; 

for  (int  j  1  =  0;  j  1  <=  n_rad;  j  1  ++)  { 


if  (jl=0){n  =  0;} 

else  if  (jl=l){n  =  1;} 

else  if  (j  1  >  l){n  =  r_fil[jl-2]+  n;} 

for  (int  m  =  0;  m  <=  on_node[n]-l;  m++){ 
for  (int  k  =  0;  k  <=  n  fil- 1 ;  k++)  { 
for  (int  j  =  0;  j  <=  on_node[k]-2;  j++){ 


double  gam  eq; 

double  dlx  =  x_n[j+l][k]-x_n[j][k]; 
double  dly  =  y_n[j+l][k]-y_n[j][k]; 


148 


double  dlz  =  z_n[j+l][k]-z_n|j][k]; 

double  rx  =  x_n[m][n]  -  (0.5*(x_n[j+l][k]  -  x_n[j][k])  +  x_n[j][k]); 
double  ry  =  y_n[m][n]  -  (0.5*(y_n[j+l][k]  -  y_n[j][k])  +  y_n[j][k]); 
double  rz  =  z_n[m][n]  -  (0.5*(z_n[j+l][k]  -  z_n[j][k])  +  z_n[j][k]); 

double  rcrossdrx  =  rz*dly  -  ry*dlz; 
double  r  cross  dr  y  =  -(rz*dlx  -  rx*dlz); 
double  r_cross_dr_z  =  ry*dlx  -  rx*dly; 

double  absrk  =  sqrt(rx*rx  +  ry*ry  +  rz*rz); 

double  s  =  l/(pow(l+(alpha*(sigma[m][n]*sigma[m][n]  +  sigma[j][k]*sigma|j][k])/ 
(2*abs_rk*abs_rk)),  1.5)); 

if  (0.5*(z_n[j+l][k]  _  z_n[j][k])  +  z_n[j][k]  <=  0){ 
gameq  =  0.0; 

} 

else  { 

gameq  =  gamma[j][k]; 

} 

double  u_ind_dummy_x  =  (-gam_eq/(4.0*pi))*r_cross_dr_x*s/(pow(abs_rk,3)); 
double  u_ind_dummy_y  =  (-gam_eq/(4.0*pi))*r_cross_dr_y*s/(pow(abs_rk,3)); 
double  u_ind_dummy_z  =  (-gam_eq/(4.0*pi))*r_cross_dr_z*s/(pow(abs_rk,3)); 

if  (j  ==  0  &&  k  ==  0)  { 

uinducedx  =  uinddummyx; 
uinducedy  =  uinddummyy; 
uinducedz  =  uinddummyz; 

} 

else  { 

uinducedx  =  uinducedx  +  u_ind_dummy_x; 
uinducedy  =  uinducedy  +  uinddummyy; 
uinducedz  =  u_induced_z  +  uinddummyz; 

} 


} 

} 

if  (z_n[m][n]  <  0){ 
u_induced_total_x[m][n]  =  0; 
u_induced_total_y[m][n]  =  0; 
u_induced_total_z[m][n]  =  0; 

} 

else  { 

u_induced_total_x[m][n]  =  u_induced_total_x[m][n]  +  uinducedx; 
u_induced_total_y[m][n]  =  u_induced_total_y[m][n]  +  uinducedy; 
u_induced_total_z[m][n]  =  u_induced_total_z[m][n]  +  uinducedz; 

}  '  '  ~ 

} 
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} 

//calculate  induced  velocities  on  and  from  radial  spoke  calculation 
r_calc_ind(); 

//calculate  induced  velocity  for  all  nodes,  based  on  symmetry 

sym_calc(); 

r_sym_calc(); 

uupdateQ; 


void  r_calc_ind() 

{ 

uinducedx  -  0; 
uinducedy  =  0; 
u  induced  z  -  0; 

r_on_fil  =  0; 

for  (int  n  =  0;  n  <=  rnfil- 1 ;  n++)  { 
if  (r_z_n[0][n]  >=  0)  { 
r_on_fil  -  ronfil  +  1  ; 

} 

} 

if  (r_on_fil  >=  r_n_fil){ 
ronfil  =  rnfil; 

} 

//calculate  velocity  induced  on  vortex  tube  filaments  by  spokes 
int  n; 

for  (int  j  1  =  0;  j  1  <=  n  rad;  j  1++)  { 


if  (jl=0){n  =  0;} 

else  if  (jl=l){n  =  1;} 

else  if  (jl  >  l){n  =  r_fil[jl-2]+  n;} 

for  (int  m  =  0;  m  <=  on_node[n]-l ;  m++) { 

for  (int  k  =  0;  k  <=  r  on  fil- 1 ;  k++)  { 
for  (int  j  =  0;  j  <=  r_max_node-2;  j++){ 


double  gameq; 

double  dlx  =  r_x_n[j+l][k]-r_x_n|j][k]; 
double  dly  =  r_y_n[j+l][k]-r_y_n[j][k]; 
double  dlz  =  r_z_n[j+ 1  ] [k]-r_z_n[j ] [k] ; 

double  rx  =  x_n[m][n]  -  (0.5*(r_x_n[j+l][k]  -  r_x_n[j][k])  +  r_x_n[j][k]); 
double  ry  =  y_n[m][n]  -  (0.5*(r_y_n[j+l][k]  -  r_y_n[j][k])  +  r_y_n[j][k]); 
double  rz  =  z_n[m][n]  -  (0.5*(r_z_n[j+l][k]  -  r_z_n[j][k])  +  r_z_n[j][k]); 


double  r  cross  dr  x  =  rz*dly  -  ry*dlz; 
double  r  cross  dr  y  =  -(rz*dlx  -  rx*dlz); 
double  r_cross_dr_z  =  ry*dlx  -  rx*dly; 
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double  absrk  =  sqrt(rx*rx  +  ry*ry  +  rz*rz); 

double  s  =  l/(pow(l+(alpha*(sigma[m][n]*sigma[m][n]  +  r_sigma[j][k]* 
r_sigma[j]  [k])/(2*abs_rk*abs_rk))?  1.5)); 

if  (0.5*(r_z_n[j+l][k]  -  r_z_n[j][k])  +  r_z_n[j][k]  <=  0){ 
gameq  =  0.0; 

} 

else  { 

gameq  =  -rgamma[j][k]; 

} 

double  uinddummyx  =  (-gam_eq/(4.0*pi))*r_cross_dr_x*s/(pow(abs_rk,3)); 
double  uinddummyy  =  (-gam_eq/(4.0*pi))*r_cross_dr_y*s/(pow(abs_rk,3)); 
double  uinddummyz  =  (-gam_eq/(4.0*pi))*r_cross_dr_z*s/(pow(abs_rk,3)); 

if  (j  =  0  &&  k  ==  0)  { 

uinducedx  =  uinddummyx; 
uinducedy  =  uinddummyy; 
uinducedz  =  u_i  n  dd  u  m  m  y_z ; 

} 

else  { 

uinducedx  =  uinducedx  +  uinddummyx; 
uinducedy  =  uinducedy  +  uinddummyy; 
uinducedz  =  u_induced_z  +  uinddummyz; 

} 


} 

} 

if  (z_n[m][n]  <  0){ 
u_induced_total_x[m][n]  =  0; 
u_induced_total_y[m][n]  =  0; 
u_induced_total_z[m][n]  =  0; 

} 

else  { 

u_induced_total_x[m][n]  =  u_induced_total_x[m][n]  +  uinducedx; 
u_induced_total_y[m][n]  =  u_induced_total_y[m][n]  +  u_induced_y; 
u_induced_total_z[m][n]  =  u_induced_total_z[m][n]  +  u_induced_z; 

}  '  "  . " 

} 

} 

//calculate  velocity  induced  on  spokes  by  vortex  tube  filaments 

for  (int  n  =  0;  n  <=  r  n  fil- 1 ;  n++)  { 
for  (int  m  =  0;  m  <=  r_max_node-l ;  m++)  { 
r_u_induced_total_x[m][n]  =  0; 
r_u_induced_total_y[m][n]  =  0; 
r_u_induced_total_z[m]  [n]  =  0; 
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} 

} 

uinducedx  -  0; 
uinducedy  =  0; 
uinducedz  =  0; 


for  (int  j  1  =  0;  j  1  <=  r_on_fil/n_spokes- 1 ;  j  1++)  { 
n  =  n_spokes*j  1 ; 

for  (int  m  =  0;  m  <=  r_max_node-l ;  m++){ 

for  (int  k  =  0;  k  <=  n  fil- 1 ;  k++)  { 
for  (int  j  =  0;  j  <=  on_node[k]-2;  j++){ 


double  gameq; 

double  dlx  =  x_n[j+l][k]-x_n[j][k]; 
double  dly  =  y_n[j+l][k]-y_n[j][k]; 
double  dlz  =  z_n[j+l][k]-z_n[j][k]; 

double  rx  =  r_x_n[m][n]  -  (0.5*(x_n[j+l][k]  -  x_n[j][k])  +  x_n[j][k]); 
double  ry  =  r_y_n[m][n]  -  (0.5*(y_n[j+l][k]  -  y_n[j][k])  +  y_n[j][k]); 
double  rz  =  r_z_n[m][n]  -  (0.5*(z_n[j+l][k]  -  z_n[j][k])  +  z_n[j][k]); 

double  rcrossdrx  =  rz*dly  -  ry*dlz; 
double  r  cross  dr  y  =  -(rz*dlx  -  rx*dlz); 
double  r  cross  dr  z  =  ry*dlx  -  rx*dly; 

double  abs  rk  =  sqrt(rx*rx  +  ry*ry  +  rz*rz); 

double  s  =  l/(pow(l+(alpha*(r_sigma[m][n]*r_sigma[m][n]  + 
sigma[j][k]*sigma|j][k])/(2*abs_rk*abs_rk)),1.5)); 

if  (0.5*(z_n[j+l][k]  -  z_n[j][k])  +  z_n[j][k]  <=  0){ 
gameq  =  0.0; 

} 

else  { 

gameq  =  gamma[j][k]; 

} 

double  u_ind_dummy_x  =  (-gam_eq/(4.0*pi))*r_cross_dr_x*s/(pow(abs_rk,3)); 
double  u_ind_dummy_y  =  (-gam_eq/(4.0*pi))*r_cross_dr_y*s/(pow(abs_rk,3)); 
double  u_ind_dummy_z  =  (-gam_eq/(4.0*pi))*r_cross_dr_z*s/(pow(abs_rk,3)); 

if  (j  ==  0  &&  k  ==  0)  { 

uinducedx  =  uinddummyx; 
uinducedy  =  uinddummyy; 
uinducedz  =  uinddummyz; 

} 

else  { 
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uinducedx  -  uinducedx  +  u_ind_dummy_x; 
uinducedy  =  uinducedy  +  uinddummyy; 
uinducedz  =  u_induced_z  +  uinddummyz; 
} 


} 

} 

if  (r_z_n[m][n]  <  0){ 
r_u_mduced_total_x[m][n]  -  0; 
r_u_mduced_total_y[m][n]  -  0; 
r_u_mduced_total_z[m][n]  -  0; 

} 

else  { 

r_u_mduced_total_x[m]  [n]  =  r_u_induced_total_x[m]  [n]+u_induced_x; 
r  ii  induced  total  y [m]  [n]  -  r_u_induced_total_y[m][n]+u_induced_y; 
r_u_induced_total_z[m]  [n]  -  r_u_induced_total_z[m]  [n]+u_induced_z; 
} 

} 

} 

} 

// - Function  init  ind - 

//  Initialize  induced  velocity  calculations 

// - 

void  init_ind() 

{ 

for  (int  n  =  0;  n  <=  n_fil- 1 ;  n++)  { 
on_node[n]  -  0; 

for  (int  m  =  0;  m  <=  max  node-l ;  m++)  { 


u_induced_total_x[m]  [n]  -  0; 
u_induced_total_y [m]  [n]  -  0; 
u_induced_total_z[m][n]  -  0; 

if  (z_n[m][n]  >=  0)  { 

on_node[n]  =  on_node[n]  +  1; 

} 


if  (on_node[n]  >=  max_node){ 
on_node[n]  -  max  node; 

} 

} 

} 

// - Function  init  ind - 

//  Update  node  positions  and  filament  parameters 

// - 

void  u_update()  { 

double  **seg_length  =  new  double  *[max_node]; 


for  (int  i  =  O;i<max_node;i++){seg_length[i]=  new  double  [n  fil];} 

for  (int  n  =  0;  n  <=  n  fil- 1 ;  n++)  { 

for  (int  m  =  0;  m  <=  max_node-l ;  m++)  { 

x_n[m][n]  -  x_n[m][n]  +  u_induced_total_x[m][n]*dt; 

y_n[m][n]  -  y_n[m][n]  +  u_induced_total_y[m][n]*dt; 

z_n[m][n]  =  z_n[m][n]  +  (u_induced_total_z[m][n]  +  U_instant[t_index])*dt; 

} 

} 

for  (int  n  =  0;  n  <=  n_fil-l;  n++){ 
for  (int  m  =  0;  m  <=  max_node-2;  m++){ 
seg_length[m][n]  =  sqrt((x_n[m+l][n]-x_n[m][n])* 

(x_n[m+ 1  ]  [n]-x_n[m]  [n])+(y_n[m+ 1  ]  [n]-y_n[m]  [n])*(y_n[m+ 1  ]  [n]- 
y_n[m]  [n])+(z_n[m+ 1  ]  [n]-z_n[m]  [n])*(z_n[m+l  ]  [n]-z_n[m]  [n])); 

sigma[m][n]  =  sqrt((core_con[n]/seg_length[m][n])); 

} 

} 

//update  radial  spokes 

r  seg-length  -  new  double  *[r_max_node]; 

for  (int  i  =  O;i<r_max_node;i++){r_seg_length[i]=  new  double  [r  n  fil] ; } 

for  (int  n  =  0;  n  <=  r  n  fil- 1 ;  n++)  { 

for  (int  m  =  0;  m  <=  r  max  node- 1 ;  m++)  { 

r_x_n[m][n]  =  r_x_n[m][n]  +  r_u_induced_total_x[m][n]*dt; 
r_y_n[m][n]  =  r_y_n[m][n]  +  r_u_induced_total_y[m][n]*dt; 
r_z_n[m][n]  =  r_z_n[m][n]  +  (r_u_induced_total_z[m][n]  +  U_instant[t_index])*dt; 
} 

} 

for  (int  n  =  0;  n  <=  r  n  fil- 1 ;  n++)  { 
for  (int  m  =  0;  m  <=  r_max_node-2;  m++){ 
r_seg_length[m][n]  -  sqrt((r_x_n[m+l][n]-r_x_n[m][n])* 

(r_x_n[ni+l  ]  [n]-r_x_n[m]  [n])+(r_y_n[m+l  ]  [n]-r_y_n[m]  [n])* 

(r_y_n[m+l  ]  [n]-r_y_n[m]  [n])+(r_z_n[m+l  ]  [n]-r_z_n[m]  [n])* 

(r_z_n  [m+ 1  ]  [n]  -r_z_n  [m]  [n] )) ; 

r_sigma[m][n]  =  sqrt((r_core_con[m][n]/r_seg_length[m][n])); 

} 

} 


//cout  «"Umax:  ”«Umax«M  Umin:  "«Umin«endl; 
delete []  seg_length; 

} 

void  sym_calc()  { 
int  n; 

for  (int  j  1  =  1;  jl  <=  n  rad;  jl++){ 
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if  01 — l)(n=  1;} 

else  if  (j  1  >  1) {n  =  r_fil[jl-2]+  n;} 

for  (int  k  =  n;  k  <=  n+r_fil[j  1  - 1  ]- 1 ;  k++)  { 
for  (int  m  =  0;  m  <=  on_node[k]-l;  m++){ 

double  urad  =  (x_n[m][n]/sqrt(x_n[m][n]*x_n[m][n]+y_n[m][n]* 
y_n[m]  [n]))*u_induced_total_x[m]  [n]+(y_n[m]  [n]/sqrt(x_n[m]  [n]  * 
x_n[m]  [n]+y_n[m]  [n]  *y_n[m]  [n]))*u_induced_total_y[m]  [n] ; 
double  u_theta  =  -(y_n[m][n]/sqrt(x_n[m][n]*x_n[m][n]+y_n[m][n]* 
y_n[m]  [n]))*u_induced_total_x[m]  [n]+(x_n[m]  [n]/sqrt(x_n[m]  [n]  * 
x_n[m]  [n]+y_n[m]  [n]  *y_n[m]  [n]))*u_induced_total_y[m]  [n] ; 
double  u_ax  =  u_induced_total_z[m][n]; 

u_induced_total_x[m]  [k]=  (x_n[m]  [k]/sqrt(x_n[m]  [k]  *x_n[m]  [k]+ 
y_n[m]  [k]  *y_n[m]  [k]))*(u_rad)-  (y_n[m]  [k]/sqrt(x_n[m]  [k]  * 
x_n[m]  [k]+y_n[m]  [k]  *y_n[m]  [k]))*(u_theta); 
u_induced_total_y  [m]  [k]=  (y_n[m]  [k]/sqrt(x_n[m]  [k]  *x_n[m]  [k]+ 
y_n[m]  [k]  *y_n[m]  [k]))*(u_rad)+(x_n[m]  [k]/sqrt(x_n[m]  [k]  * 
x_n[m]  [k]+y_n[m]  [k]  *y_n[m]  [k]))*(u_theta); 
u_induced_total_z [m]  [k] -  u_ax; 

} 

} 

} 

} 

void  r_sym_calc()  { 
int  n; 

for  (int  j  1  =  0 ;  j  1  <=  r  on  fil/n  spokes- 1 ;  j  1  ++)  { 
n  =  n_spokes*j  1 ; 

for  (int  k  =  n;  k  <=  n+n_spokes-l ;  k++){ 
for  (int  m  =  0;  m  <=  r_max_node-l ;  m++){ 

double  u  rad  -  (r_x_n[m][n]/sqrt(r_x_n[m][n]*r_x_n[m][n]+r_y_n[m][n]* 
r_y_n[m]  [n]))*r_u_induced_total_x[m]  [n]+(r_y_n[m]  [n]/ 
sqrt(r_x_n[m]  [n]  *r_x_n[m]  [n]+r_y_n[m]  [n]  *r_y_n[m]  [n]))* 
r_u_induced_total_y  [m]  [n] ; 

double  u_theta  =  -(r_y_n[m][n]/sqrt(r_x_n[m][n]*r_x_n[m][n]+ 
r_y_n[m]  [n]  *r_y_n[m]  [n]))*r_u_induced_total_x[m]  [n]+(r_x_n[m]  [n]/ 
sqrt(r_x_n[m]  [n]  *r_x_n[m]  [n]+r_y_n[m]  [n]  *r_y_n[m]  [n]))* 
ruinducedtotaly  [m]  [n] ; 
double  u_ax  =  r_u_induced_total_z[m][n]; 

r_u_induced_total_x[m]  [k]=  (r_x_n[m]  [k]/sqrt(r_x_n[m]  [k]  * 
r_x_n[m]  [k]+r_y_n[m]  [k]  *r_y_n[m]  [k]))*(u_rad)-(r_y_n[m]  [k]/ 
sqrt(r_x_n[m]  [k]  *r_x_n[m]  [k]+r_y_n[m]  [k]  *r_y_n[m]  [k]))*(u_theta); 
r_u_induced_total_y  [m]  [k]=  (ryn[m]  [k]/sqrt(r_x_n[m]  [k]  * 
r_x_n[m]  [k]+r_y_n[m]  [k]  *r_y_n[m]  [k]))*(u_rad)+(r_x_n[m]  [k]/ 
sqrt(r_x_n[m]  [k]  *r_x_n[m]  [k]+r_y_n[m]  [k]  *r_y_n[m]  [k]))*(u_theta); 
ruinducedtotalzfm]  [k]=  u_ax; 

} 

} 


} 

} 


155 


D.2  Vlinec_rad_spoke_wfj  Input  Parameter  Code 

// 

//  This  is  the  input  file  for  vlinec.exe  Do  not  modify  the  structure. 

//To  change  an  input  value,  just  replace  the  existing  value 
//  Do  not  delete  the  dollar  sign  symbols  below 

// 

$ 

1  //Simulation  initial  condition  case,  see  init_incase_rad_spoke.cpp 

2.01  //t  final.  .  .  final  simulation  time 

5  //U  inf.  .  .freestream  velocity 

0.1  //dz.  .  .  initial  axial  separation  between  vortex  filament  nodes 

0.01  //dt.  .  .simulation  time  step 

101  //maxnode.  .  maxnumber  of  nodes  in  a  filament  (excluding  radial  spokes) 

16  //numer  of  filaments  at  a  given  radius 

9  //number  of  radii  (excluding  centerline) 

0  //t_out_start.  .  .  t  index  output  starts 

1  //out  skip.  .  .  number  of  iterations  before  writing  data  to  file 

0.42  //r  tubeO.  .  .  initial  radial  location  of  an  axial  vortex  filament 

0.36  //rtubel 

0.3  //r_tube2 

0.24  //r_tube3 

0.18  //r_tube4 

0.12  //r_tube5 

0.08  //r_tube6 

0.06  //r_tube7 

0.03  //r_tube8 

0.0  //r_tube9 

0.0  //r_tubel0 

0.0  //rjubell 

0.0  //r_tubel2 

0.12  //sigma_init.  .  .  initial  vortex  filament  radius 

20.0  //Cmax.  .  .  final  value  of  total  circulation  in  vortex  tube 
-0.01  //start  C  inc.  .  .  axial  location  where  linear  C  increase  starts  (assumes  minor 

-0.51  //stop  C  inc  .  .  .  axial  location  where  linear  C  increase  ends 

(avg  vel  of  2.5  *  .2) 

0.0  //Cmin.  .  .  intial  value  of  total  circulation  in  vortex  tube 

0.413  //alpha.  .  .  parameter  defining  velocity  near  filament  center 

0.0  //z_on.  .  .  axial  location  where  self  induction  starts 

0.0  //z_init.  .  .  sets  overvall  initial  location  of  vortex  tube 

0.0  //x  init.  .  .  sets  overvall  initial  location  of  vortex  tube 

0.0  //y  init.  .  .  sets  overvall  initial  location  of  vortex  tube 


movement) 


